diff --git a/.github/workflows/black.yml b/.github/workflows/black.yml deleted file mode 100644 index 09dd075..0000000 --- a/.github/workflows/black.yml +++ /dev/null @@ -1,13 +0,0 @@ -name: Black Formatting - -on: [push, pull_request] - -jobs: - linter_name: - name: runner / black formatter - runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v2 - - uses: rickstaa/action-black@v1 - with: - black_args: ". --check" \ No newline at end of file diff --git a/.github/workflows/ruff.yml b/.github/workflows/ruff.yml index edaa44a..3707d19 100644 --- a/.github/workflows/ruff.yml +++ b/.github/workflows/ruff.yml @@ -5,6 +5,22 @@ on: [push, pull_request] jobs: ruff: runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.10"] steps: - uses: actions/checkout@v3 - - uses: chartboost/ruff-action@v1 \ No newline at end of file + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v4 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install ruff codespell tomli + - name: Run Ruff + run: | + ruff check . + - name: Spelling check with codespell + run: | + codespell --toml pyproject.toml \ No newline at end of file diff --git a/.github/workflows/yapf.yml b/.github/workflows/yapf.yml new file mode 100644 index 0000000..993b996 --- /dev/null +++ b/.github/workflows/yapf.yml @@ -0,0 +1,23 @@ +name: Yapf Format + +on: [push, pull_request] + +jobs: + yapf: + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.10"] + steps: + - uses: actions/checkout@v3 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v4 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install yapf toml + - name: Run Yapf + run: | + yapf --diff --recursive . \ No newline at end of file diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml deleted file mode 100644 index 6f3ab3a..0000000 --- a/.pre-commit-config.yaml +++ /dev/null @@ -1,22 +0,0 @@ -repos: - - repo: https://github.com/ambv/black - rev: 22.3.0 - hooks: - - id: black - - repo: https://github.com/nbQA-dev/nbQA - rev: 1.2.3 - hooks: - - id: nbqa-black - - id: nbqa-isort - - id: nbqa-flake8 - - repo: https://github.com/PyCQA/isort - rev: 5.12.0 - hooks: - - id: isort - args: ["--profile", "black", "--filter-files"] - - repo: https://github.com/astral-sh/ruff-pre-commit - # Ruff version. - rev: v0.1.6 - hooks: - # Run the linter. - - id: ruff \ No newline at end of file diff --git a/README.md b/README.md index 66a148b..7b09bea 100644 --- a/README.md +++ b/README.md @@ -7,7 +7,7 @@ Language: Python License WarpX Slack -Code style: black +Code style: yapf

@@ -26,7 +26,7 @@ Currently these "environments" are all implemented using the [Firedrake](https:/ - High-level: `hydrogym.env.FlowEnv` classes implement the OpenAI `gym.Env` interface - Intermediate: Typical CFD interface with `hydrogym.FlowConfig` and `hydrogym.TransientSolver` classes - Low-level: Access to linearized operators and sparse scipy or PETSc CSR matrices -* __Modeling and anlysis tools:__ Global stability analysis (via SLEPc) and modal decompositions (via modred) +* __Modeling and analysis tools:__ Global stability analysis (via SLEPc) and modal decompositions (via modred) * __Scalable:__ Individual environments parallelized with MPI with a **highly scalable [Ray](https://github.com/ray-project/ray) backend reinforcement learning training**. # Installation @@ -39,6 +39,9 @@ This means that the latest release of Hydrogym can be simply installed via [PyPI pip install hydrogym ``` +> BEWARE: The pip-package is currently behind the main repository, and we strongly urge users to build HydroGym +> directly from the source code. Once we've stabilized the package, we will update the pip package in turn. + However, the package assumes that the solver backend is available, so in order to run simulations locally you will need to _separately_ ensure the solver backend is installed (again, currently all the environments are implemented with Firedrake). Alternatively (and this is important for large-scale RL training), the core Hydrogym package can (or will soon be able to) launch reinforcement learning training on a Ray-cluster without an underlying Firedrake install. @@ -79,7 +82,7 @@ For more detail, check out: There are currently a number of main flow configurations, the most prominent of which are: -- Periodic cyclinder wake at Re=100 +- Periodic cylinder wake at Re=100 - Chaotic pinball at Re=130 - Open cavity at Re=7500 - Backwards-facing step at Re=600 diff --git a/docs/_templates/custom-class-template.rst b/docs/_templates/custom-class-template.rst new file mode 100644 index 0000000..d64b80d --- /dev/null +++ b/docs/_templates/custom-class-template.rst @@ -0,0 +1,34 @@ +{{ fullname | escape | underline}} + +.. currentmodule:: {{ module }} + +.. autoclass:: {{ objname }} + :members: + :show-inheritance: + :inherited-members: + :special-members: __call__, __add__, __mul__ + + {% block methods %} + {% if methods %} + .. rubric:: {{ _('Methods') }} + + .. autosummary:: + :nosignatures: + {% for item in methods %} + {%- if not item.startswith('_') %} + ~{{ name }}.{{ item }} + {%- endif -%} + {%- endfor %} + {% endif %} + {% endblock %} + + {% block attributes %} + {% if attributes %} + .. rubric:: {{ _('Attributes') }} + + .. autosummary:: + {% for item in attributes %} + ~{{ name }}.{{ item }} + {%- endfor %} + {% endif %} + {% endblock %} \ No newline at end of file diff --git a/docs/_templates/custom-module-template.rst b/docs/_templates/custom-module-template.rst new file mode 100644 index 0000000..dd90e32 --- /dev/null +++ b/docs/_templates/custom-module-template.rst @@ -0,0 +1,66 @@ +{{ fullname | escape | underline}} + +.. automodule:: {{ fullname }} + + {% block attributes %} + {% if attributes %} + .. rubric:: Module attributes + + .. autosummary:: + :toctree: + {% for item in attributes %} + {{ item }} + {%- endfor %} + {% endif %} + {% endblock %} + + {% block functions %} + {% if functions %} + .. rubric:: {{ _('Functions') }} + + .. autosummary:: + :toctree: + :nosignatures: + {% for item in functions %} + {{ item }} + {%- endfor %} + {% endif %} + {% endblock %} + + {% block classes %} + {% if classes %} + .. rubric:: {{ _('Classes') }} + + .. autosummary:: + :toctree: + :template: custom-class-template.rst + :nosignatures: + {% for item in classes %} + {{ item }} + {%- endfor %} + {% endif %} + {% endblock %} + + {% block exceptions %} + {% if exceptions %} + .. rubric:: {{ _('Exceptions') }} + + .. autosummary:: + :toctree: + {% for item in exceptions %} + {{ item }} + {%- endfor %} + {% endif %} + {% endblock %} + +{% block modules %} +{% if modules %} +.. autosummary:: + :toctree: + :template: custom-module-template.rst + :recursive: +{% for item in modules %} + {{ item }} +{%- endfor %} +{% endif %} +{% endblock %} \ No newline at end of file diff --git a/docs/api.rst b/docs/api.rst new file mode 100644 index 0000000..8a394ab --- /dev/null +++ b/docs/api.rst @@ -0,0 +1,6 @@ +.. autosummary:: + :toctree: _autosummary + :template: custom-module-template.rst + :recursive: + + hydrogym \ No newline at end of file diff --git a/docs/api_reference/api.rst b/docs/api_reference/api.rst deleted file mode 100644 index ec94338..0000000 --- a/docs/api_reference/api.rst +++ /dev/null @@ -1,7 +0,0 @@ -API -=== - -.. autosummary:: - :toctree: generated - - lumache diff --git a/docs/api_reference/index.rst b/docs/api_reference/index.rst deleted file mode 100644 index ac9d53c..0000000 --- a/docs/api_reference/index.rst +++ /dev/null @@ -1,7 +0,0 @@ -API Reference -============= - -.. toctree:: - :maxdepth: 1 - - api diff --git a/docs/conf.py b/docs/conf.py index 99097a5..124ed0b 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -1,4 +1,8 @@ # Configuration file for the Sphinx documentation builder. +import os +import sys + +sys.path.insert(0, os.path.abspath("..")) # -- Project information @@ -18,10 +22,15 @@ "sphinx.ext.autosummary", "sphinx.ext.autosectionlabel", "sphinx.ext.intersphinx", + "sphinx.ext.autodoc.typehints", "sphinx.ext.viewcode", "sphinx.ext.napoleon", "sphinx.ext.mathjax", ] +autosummary_generate = True +autoclass_content = "both" +html_show_sourcelink = False +autodoc_inherit_docstrings = True intersphinx_mapping = { "rtd": ("https://docs.readthedocs.io/en/stable/", None), @@ -31,6 +40,12 @@ intersphinx_disabled_domains = ["std"] templates_path = ["_templates"] +autodoc_mock_imports = [ + "firedrake", + "pyadjoint", + "ufl", + "mpi4py", +] # -- Options for HTML output diff --git a/docs/dev_notes/index.rst b/docs/dev_notes/index.rst index d2c32ef..b9c5340 100644 --- a/docs/dev_notes/index.rst +++ b/docs/dev_notes/index.rst @@ -1,5 +1,5 @@ -Developer Notes -=============== +Developer Documentation +======================= .. toctree:: :maxdepth: 1 @@ -7,4 +7,4 @@ Developer Notes docker_development testing distributed_backend - further_backend \ No newline at end of file + further_backend diff --git a/docs/html/.buildinfo b/docs/html/.buildinfo new file mode 100644 index 0000000..7439fdc --- /dev/null +++ b/docs/html/.buildinfo @@ -0,0 +1,4 @@ +# Sphinx build info version 1 +# This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done. +config: 399fab98a0ccd21ea9e49f55b2e83d33 +tags: 645f666f9bcd5a90fca523b33c5a78b7 diff --git a/docs/index.rst b/docs/index.rst index 17fce12..6792e97 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -12,9 +12,9 @@ Core Features * **Hierarchical:** Designed for analysis and controller design from a high-level black-box interface to low-level operator access. - * High-level: `hydrogym.env.FlowEnv` classes implement the OpenAI `gym.env` interface. - * Intermediate-level: Provides a typical CFD interface with `hydrogym.FlowConfig`, and `hydrogym.TransientSolver`. - * Low-level: Enables access to linearized operators, and sparse scipy or PETSc CSR matrices. + * High-level: `hydrogym.env.FlowEnv` classes implement the OpenAI `gym.env` interface. + * Intermediate-level: Provides a typical CFD interface with `hydrogym.FlowConfig`, and `hydrogym.TransientSolver`. + * Low-level: Enables access to linearized operators, and sparse scipy or PETSc CSR matrices. * **Modeling and Analysis Tools:** Provides global stability analysis (via SLEPc) and modal decomposition (via modred). * **Scalable:** Individual environments parallelized with MPI with a **highly scalable** `Ray `_ **backend reinforcement learning training.** @@ -26,7 +26,7 @@ Core Features .. toctree:: :hidden: - :maxdepth: 2 + :maxdepth: 3 installation quickstart @@ -35,4 +35,4 @@ Core Features integrations/index glossary dev_notes/index - api_reference/index + API Reference <_autosummary/hydrogym> diff --git a/docs/quickstart.rst b/docs/quickstart.rst index 5e944df..d495f3e 100644 --- a/docs/quickstart.rst +++ b/docs/quickstart.rst @@ -4,7 +4,7 @@ Quickstart Before getting started, follow the :ref:`installation` guide to get everything set up. Running an environment ------------- +---------------------- The `overview notebook `_ has a quick tour of some of the basic ideas of how to get Hydrogym to work in code - this is an @@ -19,6 +19,7 @@ reinforcement learning. But for now we'll just use the built-in flow environments, so we need to import `hydrogym.firedrake`: .. code-block:: python + import hydrogym.firedrake as hgym The next step is to set up the "environment" (in RL terminology), which defines the geometry, boundary @@ -26,6 +27,7 @@ conditions, actuation, transient solver, etc. of the flow configuration. We do `FlowEnv` configuration options: .. code-block:: python + env_config = { "flow": hgym.Cylinder, "flow_config": { @@ -47,12 +49,13 @@ typically expect better behavior in the first few timesteps when starting from a The `FlowEnv` provides an OpenAI Gym-compliant API for flow control, so integrating the flow is as simple as: .. code-block:: python + actuation = ctrl(obs) obs, reward, done, info = env.step(actuation) Default flow configurations: ------------- +---------------------------- Hydrogym includes several simplistic flow configurations that are benchmark problems in flow control and reduced-order modeling. This section describes a little about the benchmark configuration, control, and @@ -61,7 +64,7 @@ also has a couple of useful references for the flow, although these should by no comprehensive literature reviews. Cylinder wake -********************** +************* .. image:: _static/imgs/cylinder.png :width: 600 @@ -85,7 +88,7 @@ By default, the available measurements are the lift and drag coefficients on the Fluidic pinball -********************** +*************** .. image:: _static/imgs/pinball.png :width: 600 @@ -106,7 +109,8 @@ is chaotic, making this a much more challenging control problem. * `Maceda, et al (2021) `_ Open cavity flow -********************** +**************** + .. image:: _static/imgs/cavity.png :width: 600 @@ -123,7 +127,8 @@ present in the cylinder and pinball flows. * `Callaham, Brunton, Loiseau (2021) `_ Backwards-facing step -********************** +********************* + .. image:: _static/imgs/step.png :width: 600 @@ -138,4 +143,3 @@ analysis shows the flow is most sensitive to disturbances. * `Boujo & Gallaire (2015) `_ * `Beneddine, et al (2016) `_ * `Ducimetière, et al (2022) `_ - diff --git a/docs/requirements.txt b/docs/requirements.txt index e5ad727..6eca054 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -8,7 +8,7 @@ alabaster==0.7.12 # via sphinx babel==2.10.3 # via sphinx -certifi==2022.6.15 +certifi==2023.7.22 # via requests charset-normalizer==2.1.0 # via requests @@ -20,19 +20,19 @@ idna==3.3 # via requests imagesize==1.4.1 # via sphinx -jinja2==3.1.2 +jinja2==3.1.3 # via sphinx markupsafe==2.1.1 # via jinja2 packaging==21.3 # via sphinx -pygments==2.12.0 +pygments==2.15.0 # via sphinx pyparsing==3.0.9 # via packaging pytz==2022.1 # via babel -requests==2.28.1 +requests==2.31.0 # via sphinx snowballstemmer==2.2.0 # via sphinx @@ -55,5 +55,5 @@ sphinxcontrib-qthelp==1.0.3 # via sphinx sphinxcontrib-serializinghtml==1.1.5 # via sphinx -urllib3==1.26.9 +urllib3==1.26.18 # via requests diff --git a/docs/user_docs/ray.rst b/docs/user_docs/ray.rst index 1ea27cd..bfad8e8 100644 --- a/docs/user_docs/ray.rst +++ b/docs/user_docs/ray.rst @@ -17,7 +17,7 @@ this is heavily being worked on inside of Ray, we refer to Ray's documentation f `Ray cluster on Slurm `_. Before starting up the Python environment, we then need to perform a number of configuration -steps on all machines involed in our cluster. +steps on all machines involved in our cluster. 1. All machines that are meant to be used need to contain the same Pythone environment, libraries, and dependencies. diff --git a/examples/cylinder/pd-control.py b/examples/cylinder/pd-control.py index 7d73646..0d62324 100644 --- a/examples/cylinder/pd-control.py +++ b/examples/cylinder/pd-control.py @@ -1,10 +1,26 @@ +import os + +import matplotlib.pyplot as plt import numpy as np -import psutil +import psutil # For memory tracking import scipy.io as sio +from hydrogym.firedrake.utils.pd import PDController import hydrogym.firedrake as hgym output_dir = "output" +if not os.path.exists(output_dir): + os.makedirs(output_dir) + +mesh_resolution = "medium" + +element_type = "p1p1" +velocity_order = 1 +stabilization = "gls" + +# Make sure to run the transient simulation first to generate the restart file +restart = f"{output_dir}/checkpoint_{mesh_resolution}_{element_type}.h5" +checkpoint = f"{output_dir}/pd_{mesh_resolution}_{element_type}.h5" def compute_vort(flow): @@ -23,71 +39,73 @@ def log_postprocess(flow): # hgym.io.ParaviewCallback( # interval=10, filename=f"{output_dir}/pd-control.pvd", postprocess=compute_vort # ), + # hgym.utils.io.CheckpointCallback(interval=100, filename=checkpoint), hgym.io.LogCallback( postprocess=log_postprocess, nvals=3, - interval=10, + interval=1, print_fmt="t: {0:0.2f},\t\t CL: {1:0.3f},\t\t CD: {2:0.03f},\t\t RAM: {3:0.1f}%", - filename=None, + filename=f"{output_dir}/pd-control.dat", ), ] -flow = hgym.Cylinder( +flow = hgym.RotaryCylinder( Re=100, - mesh="coarse", - restart="../demo/checkpoint-coarse.h5", + mesh=mesh_resolution, + restart=restart, callbacks=callbacks, + velocity_order=velocity_order, ) -Tf = 1000 -dt = 1e-2 -n_steps = int(Tf // dt) - -u = np.zeros(n_steps) # Actuation history -y = np.zeros(n_steps) # Lift coefficient -dy = np.zeros(n_steps) # Derivative of lift coefficient -Kp = -4.0 # Proportional gain -Kd = 0.0 # Derivative gain +k = 2.0 # Base gain +theta = 4.0 # Phase angle +kp = k * np.cos(theta) +kd = k * np.sin(theta) + +tf = 100.0 +dt = 0.01 +n_steps = int(tf // dt) + 2 + +pd_controller = PDController( + kp, + kd, + dt, + n_steps, + filter_type="bilinear", + N=20, +) def controller(t, obs): - # return np.sin(t) - - i = int(t // dt) - - # Turn on feedback control halfway through - # if i > n_steps // 2: - # u[i] = -Kp * y[i - 1] - Kd * dy[i - 1] - - u[i] = -Kp * y[i - 1] - Kd * dy[i - 1] - - CL, CD = obs - - # Low-pass filter and estimate derivative - y[i] = y[i - 1] + (dt / flow.TAU) * (CL - y[i - 1]) - dy[i] = (y[i] - y[i - 1]) / dt - - sio.savemat(f"{output_dir}/pd-control.mat", {"y": y, "u": u}) - - return u[i] + # Turn on control halfway through + if t < tf / 2: + return 0.0 + return pd_controller(t, obs) hgym.integrate( - flow, t_span=(0, Tf), dt=dt, callbacks=callbacks, controller=controller) - -# for i in range(1, n_steps): -# # Turn on feedback control halfway through -# if i > n_steps // 2: -# u[i] = -Kp * y[i - 1] - Kd * dy[i - 1] - -# # Advance state and collect measurements -# (CL, CD), _, _, _ = env.step(u[i]) - -# # Low-pass filter and estimate derivative -# y[i] = y[i - 1] + (dt / env.flow.TAU) * (CL - y[i - 1]) -# dy[i] = (y[i] - y[i - 1]) / dt + flow, + t_span=(0, tf), + dt=dt, + callbacks=callbacks, + controller=controller, + stabilization=stabilization, +) -# hg.print( -# f"Step: {i:04d},\t\t CL: {y[i]:0.4f}, \t\tCL_dot: {dy[i]:0.4f},\t\tu: {u[i]:0.4f}" -# ) -# sio.savemat(f"{output_dir}/pd-control.mat", {"y": y, "u": u}) +# Load results data +data = np.loadtxt(f"{output_dir}/pd-control.dat") + +# Plot results +t = data[:, 0] +CL = data[:, 1] +CD = data[:, 2] + +fig, axs = plt.subplots(2, 1, figsize=(7, 4), sharex=True) +axs[0].plot(t, CL) +axs[0].set_ylabel(r"$C_L$") +axs[1].grid() +axs[1].plot(t, CD) +axs[1].set_ylabel(r"$C_D$") +axs[1].grid() +axs[1].set_xlabel("Time $t$") +plt.show() diff --git a/examples/cylinder/pd-phase-sweep.py b/examples/cylinder/pd-phase-sweep.py new file mode 100644 index 0000000..9786816 --- /dev/null +++ b/examples/cylinder/pd-phase-sweep.py @@ -0,0 +1,95 @@ +import os +import numpy as np +import psutil + +from hydrogym.firedrake.utils.pd import PDController +import hydrogym.firedrake as hgym + +output_dir = "output" +if not os.path.exists(output_dir): + os.makedirs(output_dir) + +mesh_resolution = "medium" + +element_type = "p1p1" +velocity_order = 1 +stabilization = "gls" + +restart = f"{output_dir}/checkpoint_{mesh_resolution}_{element_type}.h5" +checkpoint = f"{output_dir}/pd_{mesh_resolution}_{element_type}.h5" + + +def compute_vort(flow): + return (flow.u, flow.p, flow.vorticity()) + + +def log_postprocess(flow): + CL, CD = flow.get_observations() + mem_usage = psutil.virtual_memory().available * 100 / psutil.virtual_memory( + ).total + mem_usage = psutil.virtual_memory().percent + return CL, CD, mem_usage + + +callbacks = [ + # hgym.io.ParaviewCallback( + # interval=10, filename=f"{output_dir}/pd-control.pvd", postprocess=compute_vort + # ), + hgym.io.LogCallback( + postprocess=log_postprocess, + nvals=3, + interval=1, + print_fmt="t: {0:0.2f},\t\t CL: {1:0.3f},\t\t CD: {2:0.03f},\t\t RAM: {3:0.1f}%", + filename=f"{output_dir}/phase-sweep.dat", + ), +] + +flow = hgym.RotaryCylinder( + Re=100, + mesh=mesh_resolution, + restart=restart, + callbacks=callbacks, + velocity_order=velocity_order, +) + +omega = 1 / 5.56 # Vortex shedding frequency +k = 0.1 # Base gain +ctrl_time = 10 / omega # Time for each phasor - approx 10 vortex shedding periods +n_phase = 20 +phasors = np.linspace(0.0, 2 * np.pi, n_phase) # Phasor angles + +tf = (n_phase + 1) * ctrl_time +dt = 1e-2 +n_steps = int(tf // dt) + 2 + +pd_controller = PDController( + 0.0, + 0.0, + dt, + n_steps, + filter_type="bilinear", + N=20, +) + + +def controller(t, obs): + # Loop over phase angles for phasor control + # First interval is zero actuation + pd_controller.kp = 0.0 + pd_controller.kd = 0.0 + for j in range(n_phase): + if t > (j + 1) * ctrl_time: + pd_controller.kp = k * np.cos(phasors[j]) + pd_controller.kd = k * np.sin(phasors[j]) + + return pd_controller(t, obs) + + +hgym.integrate( + flow, + t_span=(0, tf), + dt=dt, + callbacks=callbacks, + controller=controller, + stabilization=stabilization, +) diff --git a/examples/cylinder/pressure-probes.py b/examples/cylinder/pressure-probes.py index d27aba6..1c98b54 100644 --- a/examples/cylinder/pressure-probes.py +++ b/examples/cylinder/pressure-probes.py @@ -3,7 +3,6 @@ """ import os - import matplotlib.pyplot as plt import numpy as np import psutil diff --git a/examples/cylinder/run-transient.py b/examples/cylinder/run-transient.py new file mode 100644 index 0000000..1dab814 --- /dev/null +++ b/examples/cylinder/run-transient.py @@ -0,0 +1,60 @@ +import os +import psutil +import hydrogym.firedrake as hgym + +output_dir = "output" +if not os.path.exists(output_dir): + os.makedirs(output_dir) + +pvd_out = None +restart = None +mesh_resolution = "medium" + +element_type = "p1p1" +velocity_order = 1 +stabilization = "gls" + +checkpoint = f"{output_dir}/checkpoint_{mesh_resolution}_{element_type}.h5" + +flow = hgym.RotaryCylinder( + Re=100, + restart=restart, + mesh=mesh_resolution, + velocity_order=velocity_order, +) + +# Time step +tf = 300.0 +dt = 0.01 + + +def log_postprocess(flow): + mem_usage = psutil.virtual_memory().percent + CL, CD = flow.get_observations() + return CL, CD, mem_usage + + +# Set up the callback +print_fmt = "t: {0:0.2f},\t\t CL: {1:0.3f},\t\t CD: {2:0.03f}\t\t Mem: {3:0.1f}" # This will format the output +log = hgym.utils.io.LogCallback( + postprocess=log_postprocess, + nvals=3, + interval=1, + print_fmt=print_fmt, + filename=None, +) + +interval = int(100 / dt) +callbacks = [ + log, + hgym.utils.io.CheckpointCallback(interval=interval, filename=checkpoint), +] + +hgym.print("Beginning integration") +hgym.integrate( + flow, + t_span=(0, tf), + dt=dt, + callbacks=callbacks, + stabilization=stabilization, +) diff --git a/examples/cylinder/step_input.py b/examples/cylinder/step_input.py index 96586b5..81e2471 100644 --- a/examples/cylinder/step_input.py +++ b/examples/cylinder/step_input.py @@ -1,5 +1,4 @@ """Simulate step function input at Re=40""" - import os import hydrogym.firedrake as hgym diff --git a/hydrogym/core.py b/hydrogym/core.py index 1dc5b06..62569ce 100644 --- a/hydrogym/core.py +++ b/hydrogym/core.py @@ -85,7 +85,7 @@ def set_state(self, q: StateType): """Set the current state fields Should be overridden if a different assignment - mechansim is used (e.g. `Function.assign`) + mechanism is used (e.g. `Function.assign`) Args: q (StateType): State to be assigned @@ -174,19 +174,17 @@ def set_control(self, act: ArrayLike = None): self.actuators[i].state = u def advance_time(self, dt: float, act: list[float] = None) -> list[float]: - """Update the current controls state. - - May involve integrating a dynamics model rather than - directly setting the controls state. Here, if actual - control is `u` and input is `v`, effectively - `du/dt = (1/tau)*(v - u)` + """Update the current controls state. May involve integrating + a dynamics model rather than directly setting the controls state. + Here, if actual control is `u` and input is `v`, effectively + `du/dt = (1/tau)*(v - u)` Args: - act (Iterable[ArrayLike]): Action inputs - dt (float): Time step + act (Iterable[ArrayLike]): Action inputs + dt (float): Time step Returns: - Iterable[ArrayLike]: Updated actuator state + Iterable[ArrayLike]: Updated actuator state """ if act is None: act = self.control_state diff --git a/hydrogym/firedrake/envs/cavity/flow.py b/hydrogym/firedrake/envs/cavity/flow.py index f6d6092..04c1b65 100644 --- a/hydrogym/firedrake/envs/cavity/flow.py +++ b/hydrogym/firedrake/envs/cavity/flow.py @@ -14,7 +14,7 @@ class Cavity(FlowConfig): DEFAULT_REYNOLDS = 7500 DEFAULT_MESH = "fine" DEFAULT_DT = 1e-4 - DEFAULT_STABILIZATION = "gls" + DEFAULT_STABILIZATION = "none" FUNCTIONS = ("q", "qB" ) # This flow needs a base flow to compute fluctuation KE diff --git a/hydrogym/firedrake/envs/cylinder/flow.py b/hydrogym/firedrake/envs/cylinder/flow.py index d28a4d6..fd02b10 100644 --- a/hydrogym/firedrake/envs/cylinder/flow.py +++ b/hydrogym/firedrake/envs/cylinder/flow.py @@ -134,7 +134,6 @@ def shear_force(self, q: fd.Function = None) -> float: # der of velocity wrt to the unit normal at the surface of the cylinder # equivalent to directional derivative along normal: - # https://math.libretexts.org/Courses/University_of_California_Davis/UCD_Mat_21C%3A_Multivariate_Calculus/13%3A_Partial_Derivatives/13.5%3A_Directional_Derivatives_and_Gradient_Vectors#mjx-eqn-DD2v du_dn = dot(self.epsilon(u), self.n) # Get unit tangent vector diff --git a/hydrogym/firedrake/utils/__init__.py b/hydrogym/firedrake/utils/__init__.py index 151959e..ff87239 100644 --- a/hydrogym/firedrake/utils/__init__.py +++ b/hydrogym/firedrake/utils/__init__.py @@ -2,4 +2,4 @@ from .stability import stability_analysis from .utils import get_array, is_rank_zero, print, set_from_array, white_noise -# from . import io, linalg, modeling, modred_interface +# from . import io, linalg, modeling diff --git a/hydrogym/firedrake/utils/linalg.py b/hydrogym/firedrake/utils/linalg.py index 18bc2f4..0ce8f28 100644 --- a/hydrogym/firedrake/utils/linalg.py +++ b/hydrogym/firedrake/utils/linalg.py @@ -8,11 +8,8 @@ import numpy as np import ufl from firedrake import logging -from modred import PODHandles, VectorSpaceHandles from scipy import sparse -from .modred_interface import Snapshot, vec_handle_mean - if TYPE_CHECKING: from hydrogym.firedrake import FlowConfig @@ -126,90 +123,8 @@ def inner_product(u, v): return inner_product -def pod( - flow: FlowConfig, - snapshot_handles: Iterable[Snapshot], - r: int, - mass_matrix, - decomp_indices=None, - remove_mean=True, - mean_dest="mean", - atol=1e-13, - rtol=None, - max_vecs_per_node=None, - verbosity=1, - output_dir=".", - modes_dest="pod", - eigvals_dest="eigvals.dat", - pvd_dest=None, - coeffs_dest="coeffs.dat", - field_name="q", -): - """ - Args: - ``flow``: ``FlowConfig`` for the POD (used for weighted inner product) - - ``snapshots``: List of Snapshot handles - - ``r``: Number of modes to compute - - Kwargs: - ``decomp_indices``: Indices to use in method of snapshots (defaults to ``range(r)``) - - NOTE: could actually take the "flow" dependence out if we replaced the PVD with a postprocessing callback... - """ - if fd.COMM_WORLD.size > 1: - raise NotImplementedError("Not yet supported in parallel") - - # Compute actual POD - if decomp_indices is None: - decomp_indices = range(r) - - inner_product = define_inner_product(mass_matrix) - - if remove_mean: - base_vec_handle = Snapshot(f"{output_dir}/{mean_dest}") - base_vec_handle.put(vec_handle_mean(snapshot_handles)) - - # Redefine snapshots with mean subtraction - snapshot_handles = [ - Snapshot(snap.filename, base_vec_handle=base_vec_handle) - for snap in snapshot_handles - ] - logging.log(logging.DEBUG, "Mean subtracted") - - logging.log(logging.DEBUG, "Computing POD") - POD = PODHandles( - inner_product=inner_product, - max_vecs_per_node=max_vecs_per_node, - verbosity=verbosity, - ) - POD.compute_decomp(snapshot_handles, atol=atol, rtol=rtol) - - # Vector handles for storing snapshots - mode_handles = [ - Snapshot(filename=f"{output_dir}/{modes_dest}{i}") for i in range(r) - ] - POD.compute_modes(range(r), mode_handles) - - POD.put_eigvals(f"{output_dir}/{eigvals_dest}") - - # Save for visualization - if pvd_dest is not None: - pvd = fd.File(f"{output_dir}/{pvd_dest}", "w") - for i, mode in enumerate(mode_handles): - u, p = mode.get().as_function().subfunctions - pvd.write(u, p, flow.vorticity(u)) - - # Compute temporal coefficients - coeffs = POD.compute_proj_coeffs().T # If all snapshots used for POD - np.savetxt(f"{output_dir}/{coeffs_dest}", coeffs, fmt="%0.6f", delimiter="\t") - - return coeffs, mode_handles - - -def project(basis_handles, data_handles, mass_matrix): - inner_product = define_inner_product(mass_matrix) - vec_space = VectorSpaceHandles(inner_product) - coeffs = vec_space.compute_inner_product_array(basis_handles, data_handles).T - return coeffs +# def project(basis_handles, data_handles, mass_matrix): +# inner_product = define_inner_product(mass_matrix) +# vec_space = VectorSpaceHandles(inner_product) +# coeffs = vec_space.compute_inner_product_array(basis_handles, data_handles).T +# return coeffs diff --git a/hydrogym/firedrake/utils/modred_interface.py b/hydrogym/firedrake/utils/modred_interface.py deleted file mode 100644 index 6cffbc7..0000000 --- a/hydrogym/firedrake/utils/modred_interface.py +++ /dev/null @@ -1,44 +0,0 @@ -import firedrake as fd -import modred as mr -import numpy as np - -from .utils import set_from_array - - -class Snapshot(mr.VecHandle): - """Interface between time series in numpy binaries and modred""" - - def __init__(self, filename: str, base_vec_handle=None, scale=None): - super().__init__(base_vec_handle=base_vec_handle, scale=scale) - - if filename[-4:] != ".npy": - filename += ".npy" - self.filename = filename - - def _get(self, filename=None): - # with fd.CheckpointFile(self.filename, 'r') as chk: - # q = chk.load_function(self.flow.mesh, name=self.name, idx=self.idx) - if filename is None: - filename = self.filename - return np.load(filename) - - def _put(self, q, filename=None): - if filename is None: - filename = self.filename - np.save(filename, q) - - def as_function(self, flow): - if fd.COMM_WORLD.size > 1: - raise NotImplementedError - - q = fd.Function(flow.mixed_space) - set_from_array(q, self.get()) - return q - - -def vec_handle_mean(vec_handles): - """Compute the mean of the vector handles and return as SnapshotVector""" - from functools import reduce - - data = [h.get() for h in vec_handles] - return reduce(lambda x1, x2: x1 + x2, data) * (1 / len(data)) diff --git a/hydrogym/firedrake/utils/pd.py b/hydrogym/firedrake/utils/pd.py new file mode 100644 index 0000000..d694bab --- /dev/null +++ b/hydrogym/firedrake/utils/pd.py @@ -0,0 +1,74 @@ +import numpy as np +import scipy.io as sio + + +def deriv_filter(filter_type, N, dt): + # u[i] = (b[0]*y[i] + b[1]*y[i-1] - a[1]*u[i-1]) / a[0] + if filter_type == "none": + a = [dt, 0] + b = [1, -1] + elif filter_type == "forward": + a = [1, N * dt - 1] + b = [N, -N] + elif filter_type == "backward": + a = [N * dt + 1, -1] + b = [N, -N] + elif filter_type == "bilinear": + a = [N * dt + 2, N * dt - 2] + b = [2 * N, -2 * N] + else: + raise ValueError(f"Unknown filter type: {filter_type}") + return a, b + + +class PDController: + + def __init__(self, + kp, + kd, + dt, + n_steps, + filter_type="none", + N=20, + debug_file=None): + self.kp = kp + self.kd = kd + self.N = N + self.dt = dt + + self.a, self.b = deriv_filter(filter_type, N, dt) + + self.i = 0 + self.t = np.zeros(n_steps) + self.u = np.zeros(n_steps) + self.y = np.zeros(n_steps) # Filtered lift coefficient + self.dy = np.zeros(n_steps) # Filtered derivative of lift coefficient + + self.debug_file = debug_file + + def __call__(self, t, obs): + self.i += 1 + + i, dt = self.i, self.dt + u, y, dy = self.u, self.y, self.dy + a, b, N = self.a, self.b, self.N + + u[i] = -self.kp * y[i - 1] - self.kd * dy[i - 1] + + CL, CD = obs + + # Low-pass filter and estimate derivative + y[i] = y[i - 1] + (1 / N) * (CL - y[i - 1]) + # Use the filtered measurement to avoid direct feedthrough + dy[i] = (b[0] * y[i] + b[1] * y[i - 1] - a[1] * dy[i - 1]) / a[0] + + if self.debug_file is not None and i % 100 == 0: + data = { + "y": y[:i], + "dy": dy[:i], + "u": u[:i], + "t": np.arange(0, i * dt, dt), + } + sio.savemat(self.debug_file, data) + + return u[i] diff --git a/test/README.md b/test/README.md new file mode 100644 index 0000000..a1fd687 --- /dev/null +++ b/test/README.md @@ -0,0 +1,39 @@ +# Testing of HydroGym + +## Quick-Start + +To run HydroGym's tests one best pulls the `HydroGym-Env` [docker container](https://hub.docker.com/repository/docker/lpaehler/hydrogym-env/general): + +```bash +docker pull lpaehler/hydrogym-env:stable +``` + +and then launches the VSCode Devcontainer into it. At that point one has Firedrake, and +all its dependencies pre-installed. One then needs to activate the virtualenv at the +command line with + +```bash +source /home/firedrake/firedrake/bin/activate +``` + +Install HydroGym + +```bash +pip install . +``` + +And is then set up to run the tests. + +## Running Tests + +```bash +cd test && python -m pytest test_pinball.py +``` + +or to run all tests + +```bash +python -m pytest . +``` + +> The gradient tests are currently not run, and are to be run at your own risk. diff --git a/test/test_cavity.py b/test/test_cavity.py index 51e91e4..12a6ce9 100644 --- a/test/test_cavity.py +++ b/test/test_cavity.py @@ -1,15 +1,11 @@ -import firedrake_adjoint as fda from ufl import sin +import pytest import hydrogym.firedrake as hgym -def test_import_coarse(): - hgym.Cavity(mesh="coarse") - - def test_import_medium(): - hgym.Cavity(mesh="medium") + hgym.Cavity(Re=500, mesh="medium") def test_import_fine(): @@ -17,70 +13,52 @@ def test_import_fine(): def test_steady(): - flow = hgym.Cavity(Re=50, mesh="coarse") - + flow = hgym.Cavity(Re=50, mesh="medium") solver = hgym.NewtonSolver(flow) solver.solve() def test_steady_actuation(): - flow = hgym.Cavity(Re=50, mesh="coarse") + flow = hgym.Cavity(Re=50, mesh="medium") flow.set_control(1.0) - solver = hgym.NewtonSolver(flow) solver.solve() def test_integrate(): - flow = hgym.Cavity(Re=50, mesh="coarse") + flow = hgym.Cavity(Re=50, mesh="medium") dt = 1e-4 - hgym.integrate(flow, t_span=(0, 10 * dt), dt=dt) + hgym.integrate( + flow, + t_span=(0, 2 * dt), + dt=dt, + # stabilization="gls" + ) def test_control(): - flow = hgym.Cavity(Re=50, mesh="coarse") + flow = hgym.Cavity(Re=50, mesh="medium") dt = 1e-4 - solver = hgym.IPCS(flow, dt=dt) + solver = hgym.SemiImplicitBDF(flow, dt=dt) num_steps = 10 for iter in range(num_steps): flow.get_observations() - flow = solver.step(iter, control=0.1 * sin(solver.t)) + flow = solver.step(iter, control=0.1 * sin(iter * solver.dt)) def test_env(): env_config = { "flow": hgym.Cavity, "flow_config": { - "mesh": "coarse", + "mesh": "medium", "Re": 10 }, - "solver": hgym.IPCS, + "solver": hgym.SemiImplicitBDF, } env = hgym.FlowEnv(env_config) - for _ in range(10): - y, reward, done, info = env.step(0.1 * sin(env.solver.t)) - - -def test_grad(): - flow = hgym.Cavity(Re=50, mesh="coarse") - - c = fda.AdjFloat(0.0) - flow.set_control(c) - - solver = hgym.NewtonSolver(flow) - solver.solve() - - (y,) = flow.get_observations() - - dy = fda.compute_gradient(y, fda.Control(c)) - - print(dy) - assert abs(dy) > 0 - - -if __name__ == "__main__": - test_import_medium() + for i in range(10): + y, reward, done, info = env.step(0.1 * sin(i * env.solver.dt)) diff --git a/test/test_cavity_grad.py b/test/test_cavity_grad.py new file mode 100644 index 0000000..1af4bc6 --- /dev/null +++ b/test/test_cavity_grad.py @@ -0,0 +1,22 @@ +import firedrake_adjoint as fda +from ufl import sin +import pytest + +import hydrogym.firedrake as hgym + + +def test_grad(): + flow = hgym.Cavity(Re=50, mesh="medium") + + c = fda.AdjFloat(0.0) + flow.set_control(c) + + solver = hgym.NewtonSolver(flow) + solver.solve() + + (y,) = flow.get_observations() + + dy = fda.compute_gradient(y, fda.Control(c)) + + print(dy) + assert abs(dy) > 0 diff --git a/test/test_cyl.py b/test/test_cyl.py index 3ec2116..bc6f76c 100644 --- a/test/test_cyl.py +++ b/test/test_cyl.py @@ -1,14 +1,8 @@ -import time - import firedrake as fd -import firedrake_adjoint as fda import numpy as np import hydrogym.firedrake as hgym - - -def test_import_coarse(): - hgym.Cylinder(mesh="coarse") +from hydrogym.firedrake.utils.pd import PDController def test_import_medium(): @@ -30,39 +24,25 @@ def test_steady(tol=1e-3): def test_steady_rotation(tol=1e-3): - flow = hgym.Cylinder(Re=100, mesh="medium") + Tf = 4.0 + dt = 0.1 + + flow = hgym.RotaryCylinder(Re=100, mesh="medium") flow.set_control(0.1) - solver = hgym.NewtonSolver(flow) - solver.solve() + solver = hgym.SemiImplicitBDF(flow, dt=dt) + + for iter in range(int(Tf / dt)): + solver.step(iter) # Lift/drag on cylinder CL, CD = flow.compute_forces() - assert abs(CL - 0.0594) < tol - assert abs(CD - 1.2852) < tol # Re = 100 - - -""" -def test_steady_grad(): - flow = hgym.Cylinder(Re=100, mesh="coarse") - - # First test with AdjFloat - omega = fda.AdjFloat(0.1) - - flow.set_control(omega) - - solver = hgym.NewtonSolver(flow) - solver.solve() - - J = flow.evaluate_objective() - dJ = fda.compute_gradient(J, fda.Control(omega)) - - assert abs(dJ) > 0 -""" + assert abs(CL + 0.06032) < tol + assert abs(CD - 1.49) < tol # Re = 100 def test_integrate(): - flow = hgym.Cylinder(mesh="coarse") + flow = hgym.Cylinder(mesh="medium") dt = 1e-2 hgym.integrate(flow, t_span=(0, 10 * dt), dt=dt) @@ -74,37 +54,50 @@ def feedback_ctrl(y, K=0.1): def test_control(): - flow = hgym.Cylinder(mesh="coarse") + k = 2.0 + theta = 4.0 + tf = 10.0 dt = 1e-2 - solver = hgym.IPCS(flow, dt=dt) + # Initialize the flow + flow = hgym.RotaryCylinder(Re=100, mesh="medium", velocity_order=1) - num_steps = 10 - for iter in range(num_steps): - y = flow.get_observations() - hgym.print(y) - flow = solver.step(iter, control=feedback_ctrl(y)) + # Construct the PDController + pd_controller = PDController( + k * np.cos(theta), + k * np.sin(theta), + 1e-2, + int(10.0 // 1e-2) + 2, + filter_type="bilinear", + N=2, + ) + + def controller(t, obs): + if t < tf / 2: + return 0.0 + return pd_controller(t, obs) + + hgym.integrate(flow, t_span=(0, tf), dt=dt, controller=controller) def test_env(): env_config = { "flow": hgym.Cylinder, "flow_config": { - "mesh": "coarse", + "mesh": "medium", }, - "solver": hgym.IPCS, + "solver": hgym.SemiImplicitBDF, } env = hgym.FlowEnv(env_config) u = 0.0 for _ in range(10): y, reward, done, info = env.step(u) - print(y) u = feedback_ctrl(y) def test_linearize(): - flow = hgym.Cylinder(mesh="coarse") + flow = hgym.Cylinder(mesh="medium") solver = hgym.NewtonSolver(flow) qB = solver.solve() @@ -114,141 +107,32 @@ def test_linearize(): def test_act_implicit_no_damp(): - print("") - print("No Damp") - time_start = time.time() - flow = hgym.Cylinder(mesh="coarse", actuator_integration="implicit") - dt = 1e-2 - solver = hgym.IPCS(flow, dt=dt) - - assert flow.actuators[0].integration == "implicit" + flow = hgym.Cylinder(mesh="medium", actuator_integration="implicit") + # dt = 1e-2 + solver = hgym.NewtonSolver(flow) # Since this feature is still experimental, modify actuator attributes *after*= flow.actuators[0].k = 0 - - # Apply steady torque for 0.1 seconds... should match analytical solution! - tf = 0.1 # sec - torque = 0.05 # Nm - analytical_sol = [torque * tf / flow.I_CM] - - # Run sim - num_steps = int(tf / dt) - for iter in range(num_steps): - flow = solver.step(iter, control=torque) - - print(flow.control_state[0].values(), analytical_sol) - - final_torque = fd.Constant(flow.control_state[0]).values() - assert np.isclose(final_torque, analytical_sol) - - print("finished @" + str(time.time() - time_start)) + solver.solve() def test_act_implicit_fixed_torque(): - print("") - print("Fixed Torque Convergence") - time_start = time.time() - flow = hgym.Cylinder(mesh="coarse", actuator_integration="implicit") - dt = 1e-3 - solver = hgym.IPCS(flow, dt=dt) + dt = 1e-4 + + # Define the flow + flow = hgym.Cylinder(mesh="medium", actuator_integration="implicit") - assert flow.actuators[0].integration == "implicit" + # Set up the solver + solver = hgym.SemiImplicitBDF(flow, dt=dt) # Obtain a torque value for which the system converges to a steady state angular velocity tf = 0.1 * flow.TAU omega = 1.0 - torque = omega / flow.TAU # Torque to reach steady-state value of `omega` + + # Torque to reach steady-state value of `omega` + torque = omega / flow.TAU # Run sim - num_steps = int(tf / dt) + num_steps = 10 * int(tf / dt) for iter in range(num_steps): flow = solver.step(iter, control=torque) - print(flow.control_state[0].values()) - - final_torque = fd.Constant(flow.control_state[0]).values() - assert np.isclose(final_torque, omega, atol=1e-3) - - print("finished @" + str(time.time() - time_start)) - - -def isordered(arr): - if len(arr) < 2: - return True - for i in range(len(arr) - 1): - if arr[i] < arr[i + 1]: - return False - return True - - -# sin function feeding into the controller -# TODO: This is too big for a unit test - should go in examples or somewhere like that -# def test_convergence_test_varying_torque(): -# print("") -# print("Convergence Test with Varying Torque Commands") -# time_start = time.time() -# # larger dt to compensate for longer tf -# dt_list = [1e-2, 5e-3, 2.5e-3, 1e-3] -# dt_baseline = 5e-4 - -# flow = gym.flow.Cylinder(mesh="coarse", control_method="indirect") -# solver = gym.ts.IPCS(flow, dt=dt_baseline) -# tf = 1e-1 # sec -# # torque = 50 # Nm - -# # First establish a baseline -# num_steps = int(tf / dt_baseline) -# for iter in range(num_steps): -# # let it run for 3/4 of a period of a sin wave in 0.1 second -# input = np.sin(47.12 * dt_baseline) -# flow = solver.step(iter, control=input) - -# baseline_solution = flow.get_ctrl_state()[0] - -# solutions = [] -# errors = [] - -# for dt in dt_list: -# flow = gym.flow.Cylinder(mesh="coarse", control_method="indirect") -# solver = gym.ts.IPCS(flow, dt=dt) - -# num_steps = int(tf / dt) -# for iter in range(num_steps): -# input = np.sin(47.2 * dt) -# flow = solver.step(iter, control=input) -# solutions.append(flow.get_ctrl_state()[0]) -# errors.append(np.abs(solutions[-1] - baseline_solution)) - -# # assert solutions converge to baseline solution -# assert isordered(errors) - -# print("finished @" + str(time.time() - time_start)) - -# Shear force test cases (shear force not fully implemented yet) - -# def test_shearForce0(): -# flow = gym.flow.Cylinder(Re=100, mesh="coarse") -# flow.set_control(fd.Constant(0.0)) -# flow.solve_steady() -# shear_force = flow.shear_force() - -# assert np.isclose(shear_force, 0.0, rtol=1e-3, atol=1e-3) - -# def test_shearForcePos(): -# flow = gym.flow.Cylinder(Re=100, mesh="coarse") -# flow.set_control(fd.Constant(0.1)) -# flow.solve_steady() -# shear_force = flow.shear_force() - -# assert shear_force < 0 - -# def test_shearForceNeg(): -# flow = gym.flow.Cylinder(Re=100, mesh="coarse") -# flow.set_control(fd.Constant(-0.1)) -# flow.solve_steady() -# shear_force = flow.shear_force() - -# assert shear_force > 0 - -if __name__ == "__main__": - # test_steady_grad() - test_control() diff --git a/test/test_cyl_grad.py b/test/test_cyl_grad.py new file mode 100644 index 0000000..c347401 --- /dev/null +++ b/test/test_cyl_grad.py @@ -0,0 +1,95 @@ +def isordered(arr): + if len(arr) < 2: + return True + for i in range(len(arr) - 1): + if arr[i] < arr[i + 1]: + return False + return True + + +""" +def test_steady_grad(): + flow = hgym.Cylinder(Re=100, mesh="coarse") + + # First test with AdjFloat + omega = fda.AdjFloat(0.1) + + flow.set_control(omega) + + solver = hgym.NewtonSolver(flow) + solver.solve() + + J = flow.evaluate_objective() + dJ = fda.compute_gradient(J, fda.Control(omega)) + + assert abs(dJ) > 0 +""" + +# sin function feeding into the controller +# TODO: This is too big for a unit test - should go in examples or somewhere like that +# def test_convergence_test_varying_torque(): +# print("") +# print("Convergence Test with Varying Torque Commands") +# time_start = time.time() +# # larger dt to compensate for longer tf +# dt_list = [1e-2, 5e-3, 2.5e-3, 1e-3] +# dt_baseline = 5e-4 + +# flow = gym.flow.Cylinder(mesh="coarse", control_method="indirect") +# solver = gym.ts.IPCS(flow, dt=dt_baseline) +# tf = 1e-1 # sec +# # torque = 50 # Nm + +# # First establish a baseline +# num_steps = int(tf / dt_baseline) +# for iter in range(num_steps): +# # let it run for 3/4 of a period of a sin wave in 0.1 second +# input = np.sin(47.12 * dt_baseline) +# flow = solver.step(iter, control=input) + +# baseline_solution = flow.get_ctrl_state()[0] + +# solutions = [] +# errors = [] + +# for dt in dt_list: +# flow = gym.flow.Cylinder(mesh="coarse", control_method="indirect") +# solver = gym.ts.IPCS(flow, dt=dt) + +# num_steps = int(tf / dt) +# for iter in range(num_steps): +# input = np.sin(47.2 * dt) +# flow = solver.step(iter, control=input) +# solutions.append(flow.get_ctrl_state()[0]) +# errors.append(np.abs(solutions[-1] - baseline_solution)) + +# # assert solutions converge to baseline solution +# assert isordered(errors) + +# print("finished @" + str(time.time() - time_start)) + +# Shear force test cases (shear force not fully implemented yet) + +# def test_shearForce0(): +# flow = gym.flow.Cylinder(Re=100, mesh="coarse") +# flow.set_control(fd.Constant(0.0)) +# flow.solve_steady() +# shear_force = flow.shear_force() + +# assert np.isclose(shear_force, 0.0, rtol=1e-3, atol=1e-3) + +# def test_shearForcePos(): +# flow = gym.flow.Cylinder(Re=100, mesh="coarse") +# flow.set_control(fd.Constant(0.1)) +# flow.solve_steady() +# shear_force = flow.shear_force() + +# assert shear_force < 0 + +# def test_shearForceNeg(): +# flow = gym.flow.Cylinder(Re=100, mesh="coarse") +# flow.set_control(fd.Constant(-0.1)) +# flow.solve_steady() +# shear_force = flow.shear_force() + +# assert shear_force > 0 diff --git a/test/test_pinball.py b/test/test_pinball.py index 7b0cfc0..f2b9b1c 100644 --- a/test/test_pinball.py +++ b/test/test_pinball.py @@ -1,25 +1,18 @@ -import firedrake_adjoint as fda import numpy as np -import pyadjoint - import hydrogym.firedrake as hgym -def test_import_coarse(): - hgym.Pinball(mesh="coarse") - - def test_import_fine(): hgym.Pinball(mesh="fine") def test_steady(tol=1e-2): - flow = hgym.Pinball(Re=30, mesh="coarse") + flow = hgym.Pinball(Re=30, mesh="fine") solver = hgym.NewtonSolver(flow) solver.solve() - CL_target = (0.0, 0.520, -0.517) # Slight asymmetry in mesh - CD_target = (1.4367, 1.553, 1.554) + CL_target = (0.0, 0.521, -0.521) # Slight asymmetry in mesh + CD_target = (1.451, 1.566, 1.566) CL, CD = flow.compute_forces() for i in range(len(CL)): @@ -28,59 +21,29 @@ def test_steady(tol=1e-2): def test_steady_rotation(tol=1e-2): - flow = hgym.Pinball(Re=30, mesh="coarse") + flow = hgym.Pinball(Re=30, mesh="fine") flow.set_control((0.5, 0.5, 0.5)) solver = hgym.NewtonSolver(flow) solver.solve() - CL_target = (0.2718, 0.5035, -0.6276) # Slight asymmetry in mesh - CD_target = (1.4027, 1.5166, 1.5696) + CL_target = (-0.2477, 0.356, -0.6274) + CD_target = (1.4476, 1.6887, 1.4488) CL, CD = flow.compute_forces() - print(CL) - print(CD) for i in range(len(CL)): assert abs(CL[i] - CL_target[i]) < tol assert abs(CD[i] - CD_target[i]) < tol -def test_steady_grad(): - flow = hgym.Pinball(Re=30, mesh="coarse") - n_cyl = len(flow.CYLINDER) - - # Option 1: List of AdjFloats - # omega = [fda.AdjFloat(0.5) for _ in range(n_cyl)] - # control = [fda.Control(omg) for omg in omega] - - # Option 2: List of Constants - # omega = [fd.Constant(0.5) for i in range(n_cyl)] - # control = [fda.Control(omg) for omg in omega] - - # # Option 3: Overloaded array with numpy_adjoint - omega = pyadjoint.create_overloaded_object(np.zeros(n_cyl)) - control = fda.Control(omega) - - flow.set_control(omega) - - solver = hgym.NewtonSolver(flow) - solver.solve() - - J = flow.evaluate_objective() - - # fda.compute_gradient(sum(CD), control) - dJ = fda.compute_gradient(J, control) - assert np.all(np.abs(dJ) > 0) - - def test_integrate(): - flow = hgym.Pinball(mesh="coarse") + flow = hgym.Pinball(mesh="fine") dt = 1e-2 hgym.integrate(flow, t_span=(0, 10 * dt), dt=dt) def test_control(): - flow = hgym.Pinball(mesh="coarse") + flow = hgym.Pinball(mesh="fine") dt = 1e-2 # Simple opposition control on lift @@ -90,7 +53,7 @@ def feedback_ctrl(y, K=None): CL = y[:3] return K @ CL - solver = hgym.IPCS(flow, dt=dt) + solver = hgym.SemiImplicitBDF(flow, dt=dt) num_steps = 10 for iter in range(num_steps): @@ -103,9 +66,9 @@ def test_env(): "flow": hgym.Pinball, "flow_config": { "Re": 30, - "mesh": "coarse", + "mesh": "fine", }, - "solver": hgym.IPCS, + "solver": hgym.SemiImplicitBDF, } env = hgym.FlowEnv(env_config) @@ -115,19 +78,7 @@ def feedback_ctrl(y, K=None): K = -0.1 * np.ones((3, 6)) # [Inputs x outputs] return K @ y - u = np.zeros(env.flow.ACT_DIM) + u = np.zeros(len(env.flow.CYLINDER)) for _ in range(10): y, reward, done, info = env.step(u) u = feedback_ctrl(y) - - -# def test_lti(): -# flow = gym.flow.Cylinder() -# qB = flow.solve_steady() -# A, M = flow.linearize(qB, backend='scipy') -# A_adj, M = flow.linearize(qB, adjoint=True, backend='scipy') - -if __name__ == "__main__": - # test_import_coarse() - # test_steady_rotation() - test_steady_grad() diff --git a/test/test_pinball_grad.py b/test/test_pinball_grad.py new file mode 100644 index 0000000..10ec36a --- /dev/null +++ b/test/test_pinball_grad.py @@ -0,0 +1,33 @@ +import firedrake_adjoint as fda +import numpy as np +import pyadjoint + +import hydrogym.firedrake as hgym + + +def test_steady_grad(): + flow = hgym.Pinball(Re=30, mesh="coarse") + n_cyl = len(flow.CYLINDER) + + # Option 1: List of AdjFloats + # omega = [fda.AdjFloat(0.5) for _ in range(n_cyl)] + # control = [fda.Control(omg) for omg in omega] + + # Option 2: List of Constants + # omega = [fd.Constant(0.5) for i in range(n_cyl)] + # control = [fda.Control(omg) for omg in omega] + + # # Option 3: Overloaded array with numpy_adjoint + omega = pyadjoint.create_overloaded_object(np.zeros(n_cyl)) + control = fda.Control(omega) + + flow.set_control(omega) + + solver = hgym.NewtonSolver(flow) + solver.solve() + + J = flow.evaluate_objective() + + # fda.compute_gradient(sum(CD), control) + dJ = fda.compute_gradient(J, control) + assert np.all(np.abs(dJ) > 0) diff --git a/test/test_step.py b/test/test_step.py index d40d2c6..4fde567 100644 --- a/test/test_step.py +++ b/test/test_step.py @@ -1,13 +1,8 @@ -import firedrake_adjoint as fda from ufl import sin import hydrogym.firedrake as hgym -def test_import_coarse(): - hgym.Step(mesh="coarse") - - def test_import_medium(): hgym.Step(mesh="medium") @@ -17,14 +12,14 @@ def test_import_fine(): def test_steady(): - flow = hgym.Step(Re=100, mesh="coarse") + flow = hgym.Step(Re=100, mesh="medium") solver = hgym.NewtonSolver(flow) solver.solve() def test_steady_actuation(): - flow = hgym.Step(Re=100, mesh="coarse") + flow = hgym.Step(Re=100, mesh="medium") flow.set_control(1.0) solver = hgym.NewtonSolver(flow) @@ -32,58 +27,45 @@ def test_steady_actuation(): def test_integrate(): - flow = hgym.Step(Re=100, mesh="coarse") + flow = hgym.Step(Re=100, mesh="medium") dt = 1e-3 - hgym.integrate(flow, t_span=(0, 10 * dt), dt=dt, method="IPCS") + hgym.integrate( + flow, + t_span=(0, 10 * dt), + dt=dt, + ) def test_integrate_noise(): - flow = hgym.Step(Re=100, mesh="coarse") + flow = hgym.Step(Re=100, mesh="medium") dt = 1e-3 - hgym.integrate(flow, t_span=(0, 10 * dt), dt=dt, method="IPCS", eta=1.0) + hgym.integrate(flow, t_span=(0, 10 * dt), dt=dt, eta=1.0) def test_control(): - flow = hgym.Step(Re=100, mesh="coarse") + flow = hgym.Step(Re=100, mesh="medium") dt = 1e-3 - solver = hgym.IPCS(flow, dt=dt) + solver = hgym.SemiImplicitBDF(flow, dt=dt) num_steps = 10 for iter in range(num_steps): flow.get_observations() - flow = solver.step(iter, control=0.1 * sin(solver.t)) + flow = solver.step(iter, control=0.1 * sin(iter * solver.dt)) def test_env(): env_config = { "flow": hgym.Step, "flow_config": { - "mesh": "coarse", + "mesh": "medium", "Re": 100 }, - "solver": hgym.IPCS, + "solver": hgym.SemiImplicitBDF, } env = hgym.FlowEnv(env_config) - for _ in range(10): - y, reward, done, info = env.step(0.1 * sin(env.solver.t)) - - -def test_grad(): - flow = hgym.Step(Re=100, mesh="coarse") - - c = fda.AdjFloat(0.0) - flow.set_control(c) - - solver = hgym.NewtonSolver(flow) - solver.solve() - - (y,) = flow.get_observations() - - dy = fda.compute_gradient(y, fda.Control(c)) - - print(dy) - assert abs(dy) > 0 + for i in range(10): + y, reward, done, info = env.step(0.1 * sin(i * env.solver.dt)) diff --git a/test/test_step_grad.py b/test/test_step_grad.py new file mode 100644 index 0000000..0f17938 --- /dev/null +++ b/test/test_step_grad.py @@ -0,0 +1,21 @@ +import firedrake_adjoint as fda +from ufl import sin + +import hydrogym.firedrake as hgym + + +def test_grad(): + flow = hgym.Step(Re=100, mesh="coarse") + + c = fda.AdjFloat(0.0) + flow.set_control(c) + + solver = hgym.NewtonSolver(flow) + solver.solve() + + (y,) = flow.get_observations() + + dy = fda.compute_gradient(y, fda.Control(c)) + + print(dy) + assert abs(dy) > 0 diff --git a/test/test_utils.py b/test/test_utils.py deleted file mode 100644 index e69de29..0000000