BDF/EXT time-stepper

docs/conf.py

@@ -20,7 +20,7 @@
     "sphinx.ext.intersphinx",
     "sphinx.ext.viewcode",
     "sphinx.ext.napoleon",
-    "sphinx.ext.mathjax"
+    "sphinx.ext.mathjax",
 ]
 
 intersphinx_mapping = {
@@ -44,14 +44,12 @@
 html_static_path = ["_static"]
 
 html_theme_options = {
-  'repository_url': 'https://github.com/dynamicslab/hydrogym',
-  'use_repository_button': True,  # add a 'link to repository' button
-  'use_issues_button': True,  # add an 'Open an Issue' button
-  'path_to_docs': (
-    'docs'
-  ),  # used to compute the path to launch notebooks in colab
-  'prev_next_buttons_location': None,
-  'show_navbar_depth': 1,
+    "repository_url": "https://github.com/dynamicslab/hydrogym",
+    "use_repository_button": True,  # add a 'link to repository' button
+    "use_issues_button": True,  # add an 'Open an Issue' button
+    "path_to_docs": ("docs"),  # used to compute the path to launch notebooks in colab
+    "prev_next_buttons_location": None,
+    "show_navbar_depth": 1,
 }
 
 # -- Options for EPUB output

examples/cavity/run-transient.py

@@ -0,0 +1,62 @@
+import firedrake as fd
+import numpy as np
+import psutil
+
+import hydrogym.firedrake as hgym
+
+output_dir = "output"
+pvd_out = None
+Re = 7500
+restart = f"{output_dir}/{Re}_steady.h5"
+
+flow = hgym.Cavity(Re=Re, restart=restart)
+
+# Store base flow for computing TKE
+flow.qB.assign(flow.q)
+
+# Random perturbation to base flow for initialization
+rng = fd.RandomGenerator(fd.PCG64(seed=1234))
+flow.q += rng.normal(flow.mixed_space, 0.0, 1e-2)
+
+# Time step
+Tf = 1.0
+
+method = "IPCS"
+dt = 2.5e-4
+
+# method = "BDF"
+# dt = 5e-3  # 12:13
+
+
+def log_postprocess(flow):
+    KE = 0.5 * fd.assemble(fd.inner(flow.u, flow.u) * fd.dx)
+    TKE = flow.evaluate_objective()
+    CFL = flow.max_cfl(dt)
+    mem_usage = psutil.virtual_memory().percent
+    return [CFL, KE, TKE, mem_usage]
+
+
+print_fmt = (
+    "t: {0:0.3f}\t\tCFL: {1:0.3f}\t\t KE: {2:0.3e}\t\t TKE: {3:0.3e}\t\t Mem: {4:0.1f}"
+)
+callbacks = [
+    # hgym.io.ParaviewCallback(interval=100, filename=pvd_out, postprocess=compute_vort),
+    # hgym.io.CheckpointCallback(interval=1000, filename=checkpoint),
+    hgym.io.LogCallback(
+        postprocess=log_postprocess,
+        nvals=4,
+        interval=1,
+        filename=f"{output_dir}/stats.dat",
+        print_fmt=print_fmt,
+    ),
+]
+
+
+hgym.print("Beginning integration")
+hgym.integrate(
+    flow,
+    t_span=(0, Tf),
+    dt=dt,
+    callbacks=callbacks,
+    method=method,
+)

examples/cavity/solve-steady.py

@@ -3,7 +3,7 @@
 import hydrogym.firedrake as hgym
 
 output_dir = "output"
-mesh_resolution = "coarse"
+mesh_resolution = "fine"
 Re = 7500
 
 solver_parameters = {"snes_monitor": None}

examples/demo/run-transient.py

@@ -5,13 +5,14 @@
 output_dir = "."
 pvd_out = None
 restart = None
-checkpoint = "checkpoint-coarse.h5"
+checkpoint = "checkpoint.h5"
 
-flow = hgym.Cylinder(Re=100, h5_file=restart, mesh="coarse")
+flow = hgym.Cylinder(Re=100, restart=restart, mesh="medium")
 
 # Time step
 Tf = 300
-dt = 1e-2
+Tf = 1.0
+dt = 0.1
 
 
 def log_postprocess(flow):
@@ -25,15 +26,27 @@ def log_postprocess(flow):
 log = hgym.utils.io.LogCallback(
     postprocess=log_postprocess,
     nvals=3,
-    interval=10,
+    interval=1,
     print_fmt=print_fmt,
     filename=None,
 )
 
-# callbacks = [log, hgym.utils.io.CheckpointCallback(interval=100, filename=checkpoint)]
+# callbacks = [log, hgym.utils.io.CheckpointCallback(interval=10, filename=checkpoint)]
 callbacks = [
     log,
 ]
 
+
+def controller(t, y):
+    return [1.0]
+
+
 hgym.print("Beginning integration")
-hgym.integrate(flow, t_span=(0, Tf), dt=dt, callbacks=callbacks, method="IPCS")
+hgym.integrate(
+    flow,
+    t_span=(0, Tf),
+    dt=dt,
+    callbacks=callbacks,
+    method="BDF",
+    # controller=controller,
+)

examples/pinball/run-transient.py

@@ -0,0 +1,60 @@
+import numpy as np
+import psutil
+
+import hydrogym.firedrake as hgym
+
+output_dir = "."
+pvd_out = None
+restart = None
+checkpoint = "checkpoint.h5"
+
+flow = hgym.Pinball(Re=30, restart=restart, mesh="fine")
+
+# Time step
+Tf = 1.0
+
+# # user time: 2:24 (fine), 0:26 (coarse)
+# method = "IPCS"
+# dt = 1e-2
+
+# user time: 0:50 (fine), 0:07 (coarse)
+method = "BDF"
+dt = 0.1
+
+
+def log_postprocess(flow):
+    mem_usage = psutil.virtual_memory().percent
+    obs = flow.get_observations()
+    CL = obs[:3]
+    return *CL, mem_usage
+
+
+# Set up the callback
+print_fmt = "t: {0:0.2f},\t\t CL[0]: {1:0.3f},\t\t CL[1]: {2:0.03f},\t\t CL[2]: {3:0.03f}\t\t Mem: {4:0.1f}"  # noqa: E501
+log = hgym.utils.io.LogCallback(
+    postprocess=log_postprocess,
+    nvals=4,
+    interval=1,
+    print_fmt=print_fmt,
+    filename="coeffs.dat",
+)
+
+# callbacks = [log, hgym.utils.io.CheckpointCallback(interval=100, filename=checkpoint)]
+callbacks = [
+    log,
+]
+
+
+def controller(t, obs):
+    return np.array([0.0, 1.0, 1.0])
+
+
+hgym.print("Beginning integration")
+hgym.integrate(
+    flow,
+    t_span=(0, Tf),
+    dt=dt,
+    # controller=controller,
+    callbacks=callbacks,
+    method=method,
+)

examples/step/run-transient.py

@@ -0,0 +1,58 @@
+import firedrake as fd
+import numpy as np
+import psutil
+
+import hydrogym.firedrake as hgym
+
+Re = 600
+output_dir = "output"
+mesh_resolution = "fine"
+restart = f"{output_dir}/{Re}_steady.h5"
+
+flow = hgym.Step(Re=Re, mesh=mesh_resolution, restart=restart)
+
+# Store base flow for computing TKE
+flow.qB.assign(flow.q)
+
+tf = 1.0
+
+# method = "IPCS"
+# dt = 1e-3
+
+method = "BDF"
+dt = 1e-2
+
+
+def log_postprocess(flow):
+    KE = 0.5 * fd.assemble(fd.inner(flow.u, flow.u) * fd.dx)
+    TKE = flow.evaluate_objective()
+    CFL = flow.max_cfl(dt)
+    mem_usage = psutil.virtual_memory().percent
+    return [CFL, KE, TKE, mem_usage]
+
+
+print_fmt = (
+    "t: {0:0.3f}\t\tCFL: {1:0.3f}\t\t KE: {2:0.6e}\t\t TKE: {3:0.6e}\t\t Mem: {4:0.1f}"
+)
+interval = int(1e-2 / dt)
+callbacks = [
+    # hgym.io.ParaviewCallback(interval=100, filename=pvd_out, postprocess=compute_vort),
+    # hgym.io.CheckpointCallback(interval=100, filename=checkpoint),
+    hgym.io.LogCallback(
+        postprocess=log_postprocess,
+        nvals=4,
+        interval=interval,
+        filename=f"{output_dir}/{method}_stats.dat",
+        print_fmt=print_fmt,
+    ),
+]
+
+
+hgym.integrate(
+    flow,
+    t_span=(0, tf),
+    dt=dt,
+    callbacks=callbacks,
+    method=method,
+    eta=1.0,
+)

examples/step/solve-steady.py

@@ -6,6 +6,7 @@
 output_dir = "output"
 mesh_resolution = "fine"
 Re = 600
+checkpoint_prefix = f"{output_dir}/{Re}_steady"
 
 solver_parameters = {"snes_monitor": None}
 
@@ -23,7 +24,7 @@
     solver = hgym.NewtonSolver(flow, solver_parameters=solver_parameters)
     qB = solver.solve()
 
-flow.save_checkpoint(f"{output_dir}/{Re}_steady.h5")
+flow.save_checkpoint(f"{checkpoint_prefix}.h5")
 vort = flow.vorticity()
-pvd = fd.File(f"{output_dir}/{Re}_steady.pvd")
+pvd = fd.File(f"{checkpoint_prefix}.pvd")
 pvd.write(flow.u, flow.p, vort)

hydrogym/core.py

@@ -6,13 +6,15 @@
 
 class ActuatorBase:
     def __init__(self, state=0.0, **kwargs):
-        self.u = state
+        self.x = state
 
-    def set_state(self, u: float):
-        self.u = u
+    @property
+    def state(self) -> float:
+        return self.x
 
-    def get_state(self) -> float:
-        return self.u
+    @state.setter
+    def state(self, u: float):
+        self.x = u
 
     def step(self, u: float, dt: float):
         """Update the state of the actuator"""
@@ -153,14 +155,14 @@ def enlist(self, x: Any) -> Iterable[Any]:
 
     @property
     def control_state(self) -> Iterable[ActType]:
-        return [a.get_state() for a in self.actuators]
+        return [a.state for a in self.actuators]
 
     def set_control(self, act: ActType = None):
         """Directly set the control state"""
         if act is None:
             act = np.zeros(self.ACT_DIM)
         for i, u in enumerate(self.enlist(act)):
-            self.actuators[i].set_state(u)
+            self.actuators[i].state = u
 
     def update_actuators(self, act: Iterable[ActType], dt: float) -> Iterable[ActType]:
         """Update the current controls state.
@@ -176,8 +178,6 @@ def update_actuators(self, act: Iterable[ActType], dt: float) -> Iterable[ActTyp
 
         Returns:
             Iterable[ActType]: Updated actuator state
-
-        TODO: Rewrite with ActuatorBase
         """
         act = self.enlist(act)
         assert len(act) == self.ACT_DIM

hydrogym/firedrake/__init__.py

@@ -1,8 +1,8 @@
 from hydrogym.core import FlowEnv
 
 from .actuator import DampedActuator
-from .flow import FlowConfig
-from .solver import IPCS, NewtonSolver, integrate
+from .flow import FlowConfig, ScaledDirichletBC
+from .solvers import IPCS, NewtonSolver, SemiImplicitBDF, integrate
 from .utils import io, is_rank_zero, linalg, modeling, print
 
 from .envs import Cylinder, Pinball, Cavity, Step  # isort:skip