diff --git a/README.md b/README.md index 66a148b..052cc1d 100644 --- a/README.md +++ b/README.md @@ -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. diff --git a/examples/cylinder/pd-control.py b/examples/cylinder/pd-control.py index cd13ebf..e9c0427 100644 --- a/examples/cylinder/pd-control.py +++ b/examples/cylinder/pd-control.py @@ -4,8 +4,8 @@ import numpy as np import psutil # For memory tracking import scipy.io as sio -from pd import PDController +from hydrogym.firedrake.utils.pd import PDController import hydrogym.firedrake as hgym output_dir = "output" diff --git a/examples/cylinder/pd-phase-sweep.py b/examples/cylinder/pd-phase-sweep.py index 6b1fb40..46035cb 100644 --- a/examples/cylinder/pd-phase-sweep.py +++ b/examples/cylinder/pd-phase-sweep.py @@ -1,9 +1,8 @@ import os - import numpy as np import psutil -from pd import PDController +from hydrogym.firedrake.utils.pd import PDController import hydrogym.firedrake as hgym output_dir = "output" diff --git a/examples/cylinder/pressure-probes.py b/examples/cylinder/pressure-probes.py index 239ba33..aaeb3a5 100644 --- a/examples/cylinder/pressure-probes.py +++ b/examples/cylinder/pressure-probes.py @@ -2,7 +2,6 @@ Simulate the flow around the cylinder with evenly spaced surface pressure probes. """ 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 index 9b9753c..d0e2f0d 100644 --- a/examples/cylinder/run-transient.py +++ b/examples/cylinder/run-transient.py @@ -1,7 +1,5 @@ import os - import psutil - import hydrogym.firedrake as hgym output_dir = "output" diff --git a/examples/cylinder/step_input.py b/examples/cylinder/step_input.py index 98894aa..85e1923 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/firedrake/envs/cavity/flow.py b/hydrogym/firedrake/envs/cavity/flow.py index aa95cad..efc560a 100644 --- a/hydrogym/firedrake/envs/cavity/flow.py +++ b/hydrogym/firedrake/envs/cavity/flow.py @@ -1,7 +1,10 @@ import os import firedrake as fd +import matplotlib.pyplot as plt +import numpy as np import ufl +from firedrake.pyplot import tricontourf from ufl import dot, ds, grad from hydrogym.firedrake import FlowConfig, ObservationFunction, ScaledDirichletBC @@ -11,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 @@ -108,3 +111,23 @@ def evaluate_objective(self, q=None, qB=None): uB = qB.subfunctions[0] KE = 0.5 * fd.assemble(fd.inner(u - uB, u - uB) * fd.dx) return KE + + # TODO: Rendering function needs to be revisited as this is only a hot-fix + def render(self, mode="human", clim=None, levels=None, cmap="RdBu", **kwargs): + if clim is None: + clim = (-2, 2) + if levels is None: + levels = np.linspace(*clim, 10) + vort = fd.project(fd.curl(self.u), self.pressure_space) + im = tricontourf( + vort, + cmap=cmap, + levels=levels, + vmin=clim[0], + vmax=clim[1], + extend="both", + **kwargs, + ) + + cyl = plt.Circle((0, 0), 0.5, edgecolor="k", facecolor="gray") + im.axes.add_artist(cyl) diff --git a/hydrogym/firedrake/envs/cylinder/flow.py b/hydrogym/firedrake/envs/cylinder/flow.py index 5faba2d..4aa4ddb 100644 --- a/hydrogym/firedrake/envs/cylinder/flow.py +++ b/hydrogym/firedrake/envs/cylinder/flow.py @@ -131,7 +131,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/envs/pinball/flow.py b/hydrogym/firedrake/envs/pinball/flow.py index bf896a3..f04e83e 100644 --- a/hydrogym/firedrake/envs/pinball/flow.py +++ b/hydrogym/firedrake/envs/pinball/flow.py @@ -113,6 +113,7 @@ def evaluate_objective(self, q=None): CL, CD = self.compute_forces(q=q) return sum(CD) + # TODO: Needs to be revisited as the self calls here look hella suss def render(self, mode="human", clim=None, levels=None, cmap="RdBu", **kwargs): if clim is None: clim = (-2, 2) diff --git a/hydrogym/firedrake/envs/step/flow.py b/hydrogym/firedrake/envs/step/flow.py index a56fd2c..9173275 100644 --- a/hydrogym/firedrake/envs/step/flow.py +++ b/hydrogym/firedrake/envs/step/flow.py @@ -1,9 +1,11 @@ import os import firedrake as fd +import matplotlib.pyplot as plt import numpy as np import ufl from firedrake.petsc import PETSc +from firedrake.pyplot import tricontourf from ufl import dot, ds, exp, grad from hydrogym.firedrake import FlowConfig, ObservationFunction, ScaledDirichletBC @@ -143,3 +145,23 @@ def evaluate_objective(self, q=None, qB=None): uB = qB.subfunctions[0] KE = 0.5 * fd.assemble(fd.inner(u - uB, u - uB) * fd.dx) return KE + + # TODO: Rendering function needs to be revisited as this is only a hot-fix + def render(self, mode="human", clim=None, levels=None, cmap="RdBu", **kwargs): + if clim is None: + clim = (-2, 2) + if levels is None: + levels = np.linspace(*clim, 10) + vort = fd.project(fd.curl(self.u), self.pressure_space) + im = tricontourf( + vort, + cmap=cmap, + levels=levels, + vmin=clim[0], + vmax=clim[1], + extend="both", + **kwargs, + ) + + cyl = plt.Circle((0, 0), 0.5, edgecolor="k", facecolor="gray") + im.axes.add_artist(cyl) diff --git a/hydrogym/firedrake/solvers/bdf_ext.py b/hydrogym/firedrake/solvers/bdf_ext.py index dcb6259..22cc3fb 100644 --- a/hydrogym/firedrake/solvers/bdf_ext.py +++ b/hydrogym/firedrake/solvers/bdf_ext.py @@ -119,11 +119,9 @@ def _make_order_k_solver(self, k): "pc_fieldsplit_type": "schur", "pc_fieldsplit_schur_fact_type": "full", "pc_fieldsplit_schur_precondition": "selfp", - # # Default preconditioner for inv(A) # (ilu in serial, bjacobi in parallel) "fieldsplit_0_ksp_type": "preonly", - # # Single multigrid cycle preconditioner for inv(S) "fieldsplit_1_ksp_type": "preonly", "fieldsplit_1_pc_type": "hypre", diff --git a/examples/cylinder/pd.py b/hydrogym/firedrake/utils/pd.py similarity index 100% rename from examples/cylinder/pd.py rename to hydrogym/firedrake/utils/pd.py diff --git a/pyproject.toml b/pyproject.toml index c83a38e..656be1e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -71,8 +71,8 @@ black = "^23.11.0" [tool.ruff] line-length = 120 -select = ["E", "F"] -ignore = ["F401"] +lint.select = ["E", "F"] +lint.ignore = ["F401"] [build-system] requires = ["poetry-core>=1.0.0"] 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 1c18e67..a349482 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,67 +13,49 @@ 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", "Re": 10}, - "solver": hgym.IPCS, + "flow_config": {"mesh": "medium", "Re": 10}, + "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..4ee14f9 --- /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 628e71e..3dcfc68 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,144 +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 - assert flow.actuators[0].integration == "implicit" + # Define the flow + flow = hgym.Cylinder(mesh="medium", actuator_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..2e286e2 --- /dev/null +++ b/test/test_cyl_grad.py @@ -0,0 +1,98 @@ +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 17f18d7..7f739a9 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..d0bbb44 --- /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 fb9dc77..8a8f4de 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,55 +27,42 @@ 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", "Re": 100}, - "solver": hgym.IPCS, + "flow_config": {"mesh": "medium", "Re": 100}, + "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..6f577ca --- /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