diff --git a/examples/3d_examples/RestartExample/restart_example.py b/examples/3d_examples/RestartExample/restart_example.py new file mode 100644 index 0000000..0189ea2 --- /dev/null +++ b/examples/3d_examples/RestartExample/restart_example.py @@ -0,0 +1,424 @@ +import elastica as ea +import numpy as np +import sopht.utils as spu +import sopht.simulator as sps +import click +from matplotlib import pyplot as plt + + +def flow_past_rod_case( + non_dim_final_time: float, + n_elem: int, + grid_size: tuple[int, int, int], + surface_grid_density_for_largest_element: int, + cauchy_number: float, + mass_ratio: float, + froude_number: float, + stretch_bending_ratio: float, + poisson_ratio: float = 0.5, + reynolds: float = 100.0, + coupling_stiffness: float = -2e5, + coupling_damping: float = -1e2, + rod_start_incline_angle: float = 0.0, + num_threads: int = 4, + precision: str = "single", + save_data: bool = False, +) -> None: + # =================COMMON SIMULATOR STUFF======================= + grid_dim = 3 + grid_size_z, grid_size_y, grid_size_x = grid_size + real_t = spu.get_real_t(precision) + x_axis_idx = spu.VectorField.x_axis_idx() + z_axis_idx = spu.VectorField.z_axis_idx() + rho_f = 1.0 + U_free_stream = 1.0 + base_length = 1.0 + x_range = 1.8 * base_length + y_range = grid_size_y / grid_size_x * x_range + z_range = grid_size_z / grid_size_x * x_range + velocity_free_stream = [U_free_stream, 0.0, 0.0] + # =================PYELASTICA STUFF BEGIN===================== + + class RestartExampleSimulator( + ea.BaseSystemCollection, ea.Constraints, ea.Forcing, ea.Damping + ): + ... + + restart_example_simulator = RestartExampleSimulator() + start = np.array([0.2 * x_range, 0.5 * y_range, 0.75 * z_range]) + direction = np.array( + [np.sin(rod_start_incline_angle), 0.0, -np.cos(rod_start_incline_angle)] + ) + normal = np.array([0.0, 1.0, 0.0]) + base_diameter = y_range / 5.0 + base_radius = base_diameter / 2.0 + base_area = np.pi * base_radius**2 + # mass_ratio = rho_s / rho_f + rho_s = mass_ratio * rho_f + moment_of_inertia = np.pi / 4 * base_radius**4 + # Cau = (rho_f U^2 L^3 D) / EI + youngs_modulus = (rho_f * U_free_stream**2 * base_length**3 * base_diameter) / ( + cauchy_number * moment_of_inertia + ) + # Froude = g L / U^2 + gravitational_acc = froude_number * U_free_stream**2 / base_diameter + # Stretch to Bending ratio EAL2 / EI + Es_Eb = stretch_bending_ratio * moment_of_inertia / (base_area * base_length**2) + + flow_past_rod = ea.CosseratRod.straight_rod( + n_elem, + start, + direction, + normal, + base_length, + base_radius, + rho_s, + youngs_modulus=youngs_modulus, + shear_modulus=youngs_modulus / (poisson_ratio + 1.0), + ) + flow_past_rod.shear_matrix[2, 2, :] *= Es_Eb + restart_example_simulator.append(flow_past_rod) + restart_example_simulator.constrain(flow_past_rod).using( + ea.OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,) + ) + # Add gravitational forces + restart_example_simulator.add_forcing_to(flow_past_rod).using( + ea.GravityForces, acc_gravity=np.array([0.0, 0.0, -gravitational_acc]) + ) + # add damping + dl = base_length / n_elem + rod_dt = 0.01 * dl + damping_constant = 1e-3 + restart_example_simulator.dampen(flow_past_rod).using( + ea.AnalyticalLinearDamper, + damping_constant=damping_constant, + time_step=rod_dt, + ) + # =================PYELASTICA STUFF END===================== + + # ==================FLOW SETUP START========================= + kinematic_viscosity = U_free_stream * base_diameter / reynolds + flow_sim = sps.UnboundedNavierStokesFlowSimulator3D( + grid_size=grid_size, + x_range=x_range, + kinematic_viscosity=kinematic_viscosity, + with_forcing=True, + with_free_stream_flow=True, + real_t=real_t, + num_threads=num_threads, + filter_vorticity=True, + filter_setting_dict={"order": 1, "type": "multiplicative"}, + ) + # ==================FLOW SETUP END========================= + + # ==================FLOW-ROD COMMUNICATOR SETUP START====== + cosserat_rod_flow_interactor = sps.CosseratRodFlowInteraction( + cosserat_rod=flow_past_rod, + eul_grid_forcing_field=flow_sim.eul_grid_forcing_field, + eul_grid_velocity_field=flow_sim.velocity_field, + virtual_boundary_stiffness_coeff=coupling_stiffness, + virtual_boundary_damping_coeff=coupling_damping, + dx=flow_sim.dx, + grid_dim=grid_dim, + real_t=real_t, + num_threads=num_threads, + forcing_grid_cls=sps.CosseratRodSurfaceForcingGrid, + surface_grid_density_for_largest_element=surface_grid_density_for_largest_element, + ) + restart_example_simulator.add_forcing_to(flow_past_rod).using( + sps.FlowForces, + cosserat_rod_flow_interactor, + ) + # ==================FLOW-ROD COMMUNICATOR SETUP END====== + # =================TIMESTEPPING==================== + restart_example_simulator.finalize() + timestepper = ea.PositionVerlet() + do_step, stages_and_updates = ea.extend_stepper_interface( + timestepper, restart_example_simulator + ) + + if save_data: + # setup flow IO + io = spu.EulerianFieldIO( + position_field=flow_sim.position_field, + eulerian_fields_dict={ + "vorticity": flow_sim.vorticity_field, + "velocity": flow_sim.velocity_field, + }, + ) + # Initialize sphere IO + rod_io = spu.CosseratRodIO( + cosserat_rod=flow_past_rod, dim=grid_dim, real_dtype=real_t + ) + # Initialize forcing io + forcing_io = spu.IO(dim=grid_dim, real_dtype=real_t) + # Add vector field on lagrangian grid + forcing_io.add_as_lagrangian_fields_for_io( + lagrangian_grid=cosserat_rod_flow_interactor.forcing_grid.position_field, + lagrangian_grid_name="cosseratrod", + lagrangian_grid_forcing_field=cosserat_rod_flow_interactor.lag_grid_forcing_field, + lagrangiang_grid_position_mismatch_field=cosserat_rod_flow_interactor.lag_grid_position_mismatch_field, + ) + + restart_dir = "restart_data" + + spu.restart_simulation( + flow_sim=flow_sim, + restart_simulator=restart_example_simulator, + io=io, + rod_io=rod_io, + forcing_io=forcing_io, + restart_dir=restart_dir, + ) + + foto_timer = 0.0 + timescale = base_length / U_free_stream + final_time = non_dim_final_time * timescale + foto_timer_limit = 1 / 30 + time_history = [] + rod_angle = [] + data_timer = 0.0 + data_timer_limit = 0.005 * timescale + tip_time = [] + tip_position = [] + force_history = [] + + # create fig for plotting flow fields + fig, ax = spu.create_figure_and_axes() + + def rod_incline_angle_with_horizon(rod: ea.CosseratRod): + return np.rad2deg( + np.fabs( + np.arctan( + ( + rod.position_collection[z_axis_idx, -1] + - rod.position_collection[z_axis_idx, 0] + ) + / ( + rod.position_collection[x_axis_idx, -1] + - rod.position_collection[x_axis_idx, 0] + ) + ) + ) + ) + + # iterate + while flow_sim.time < final_time: + # Save data + if foto_timer > foto_timer_limit or foto_timer == 0: + foto_timer = 0.0 + ax.set_title(f"Velocity magnitude, time: {flow_sim.time / timescale:.2f}") + contourf_obj = ax.contourf( + flow_sim.position_field[x_axis_idx, :, grid_size_y // 2, :], + flow_sim.position_field[z_axis_idx, :, grid_size_y // 2, :], + np.linalg.norm( + np.mean( + flow_sim.velocity_field[ + :, :, grid_size_y // 2 - 1 : grid_size_y // 2 + 1, : + ], + axis=2, + ), + axis=0, + ), + levels=50, + extend="both", + cmap=spu.get_lab_cmap(), + ) + cbar = fig.colorbar(mappable=contourf_obj, ax=ax) + ax.scatter( + cosserat_rod_flow_interactor.forcing_grid.position_field[x_axis_idx], + cosserat_rod_flow_interactor.forcing_grid.position_field[z_axis_idx], + s=5, + color="k", + ) + spu.save_and_clear_fig( + fig, + ax, + cbar, + file_name="snap_" + str("%0.5d" % (flow_sim.time * 100)) + ".png", + ) + time_history.append(flow_sim.time) + rod_angle.append(rod_incline_angle_with_horizon(flow_past_rod)) + forces = np.sum(cosserat_rod_flow_interactor.lag_grid_forcing_field, axis=1) + force_history.append(forces.copy()) + print( + f"time: {flow_sim.time:.2f} ({(flow_sim.time/final_time*100):2.1f}%), " + f"max_vort: {np.amax(flow_sim.vorticity_field):.4f}, " + f"rod angle: {rod_incline_angle_with_horizon(flow_past_rod):2.2f}, " + f"vort divg. L2 norm: {flow_sim.get_vorticity_divergence_l2_norm():.4f}," + " grid deviation L2 error: " + f"{cosserat_rod_flow_interactor.get_grid_deviation_error_l2_norm():.6f}," + f" total force: {forces}," + f" force norm: {np.linalg.norm(forces)}" + ) + if save_data: + io.save( + h5_file_name="sopht_" + + str("%0.4d" % (flow_sim.time * 100)) + + ".h5", + time=flow_sim.time, + ) + rod_io.save( + h5_file_name="rod_" + str("%0.4d" % (flow_sim.time * 100)) + ".h5", + time=flow_sim.time, + ) + forcing_io.save( + h5_file_name="forcing_grid_" + + str("%0.4d" % (flow_sim.time * 100)) + + ".h5", + time=flow_sim.time, + ) + ea.save_state(restart_example_simulator, restart_dir, flow_sim.time) + + # save diagnostic data + if data_timer >= data_timer_limit or data_timer == 0: + data_timer = 0.0 + tip_time.append(flow_sim.time / timescale) + tip_position.append((flow_past_rod.position_collection[(x_axis_idx), -1])) + + # compute timestep + flow_dt = flow_sim.compute_stable_timestep(dt_prefac=0.5) + # flow_dt = rod_dt + + # timestep the rod, through the flow timestep + rod_time_steps = int(flow_dt / min(flow_dt, rod_dt)) + local_rod_dt = flow_dt / rod_time_steps + rod_time = flow_sim.time + for i in range(rod_time_steps): + rod_time = do_step( + timestepper, + stages_and_updates, + restart_example_simulator, + rod_time, + local_rod_dt, + ) + # timestep the cosserat_rod_flow_interactor + cosserat_rod_flow_interactor.time_step(dt=local_rod_dt) + # evaluate feedback/interaction between flow and rod + cosserat_rod_flow_interactor() + + flow_sim.time_step( + dt=flow_dt, + free_stream_velocity=velocity_free_stream, + ) + + # update timer + foto_timer += flow_dt + + # compile video + spu.make_video_from_image_series( + video_name="flow", image_series_name="snap", frame_rate=30 + ) + + # Save data + np.savetxt( + "rod_angle_vs_time.csv", + np.c_[np.array(time_history), np.array(rod_angle)], + delimiter=",", + header="time, rod_angle", + ) + + np.savetxt( + "rod_forces_vs_time.csv", + np.c_[ + np.array(time_history), + np.array(force_history), + np.linalg.norm(np.array(force_history), axis=1), + ], + delimiter=",", + header="time, force x, force y, force z, force norm", + ) + + combined_array = np.column_stack((np.array(tip_time), np.array(tip_position))) + np.savetxt("tip_displacement.csv", combined_array, delimiter=",") + + plt.figure() + plt.plot(np.array(tip_time), np.array(tip_position), label="X") + plt.xlabel("Non-dimensional time") + plt.ylabel("Tip deflection") + plt.savefig("tip_position_vs_time.png") + + +if __name__ == "__main__": + + @click.command() + @click.option("--num_threads", default=4, help="Number of threads for parallelism.") + @click.option("--nx", default=128, help="Number of grid points in x direction.") + @click.option("--u_free_stream", default=1.1, help="Free stream flow velocity.") + def simulate_flow_past_rod_with_restart( + num_threads: int, nx: int, u_free_stream: float + ) -> None: + ny = nx // 4 + nz = nx + # in order Z, Y, X + grid_size = (nz, ny, nx) + surface_grid_density_for_largest_element = nx // 8 + n_elem = 5 * nx // 16 + + click.echo(f"Number of threads for parallelism: {num_threads, }") + click.echo(f"Grid size: {nz, ny, nx ,} ") + click.echo( + f"num forcing points around the surface: {surface_grid_density_for_largest_element}" + ) + click.echo(f"num rod elements: {n_elem}") + click.echo(f"Free stream flow velocity: {u_free_stream}") + + final_time = 4 + + exp_rho_s = 1e3 # kg/m3 + exp_rho_f = 1.21 # kg/m3 + exp_youngs_modulus = 2.25e5 # Pa + exp_poisson_ratio = 0.01 + exp_base_length = 25e-3 # m + exp_base_diameter = 0.4e-3 # m + exp_kinematic_viscosity = 1.51e-5 # m2/s + exp_U_free_stream = u_free_stream # m/s + exp_gravitational_acc = 9.80665 # m/s2 + + exp_mass_ratio = exp_rho_s / exp_rho_f + exp_base_radius = exp_base_diameter / 2 + exp_base_area = np.pi * exp_base_radius**2 + exp_moment_of_inertia = np.pi / 4 * exp_base_radius**4 + exp_bending_rigidity = exp_youngs_modulus * exp_moment_of_inertia + exp_cauchy_number = ( + exp_rho_f + * exp_U_free_stream**2 + * exp_base_length**3 + * exp_base_diameter + / exp_bending_rigidity + ) + # Froude = g D / U^2 + exp_froude_number = ( + exp_gravitational_acc * exp_base_diameter / exp_U_free_stream**2 + ) + exp_Re = exp_U_free_stream * exp_base_diameter / exp_kinematic_viscosity + + # stretch to bending ratio EAL2 / EI + exp_Ks_Kb = (exp_youngs_modulus * exp_base_area * exp_base_length**2) / ( + exp_youngs_modulus * exp_moment_of_inertia + ) + + # Final deflection angle Silvaleon 2018 Eq 15 + rod_start_incline_angle = np.deg2rad(0) + + print( + "Re: {exp_Re}, Ca: {exp_cauchy_number}, Fr: {exp_froude_number}, Angle: {rod_start_incline_angle}" + ) + + flow_past_rod_case( + non_dim_final_time=final_time, + cauchy_number=exp_cauchy_number, + mass_ratio=exp_mass_ratio, + froude_number=exp_froude_number, + poisson_ratio=exp_poisson_ratio, + reynolds=exp_Re, + stretch_bending_ratio=exp_Ks_Kb, + grid_size=grid_size, + surface_grid_density_for_largest_element=surface_grid_density_for_largest_element, + n_elem=n_elem, + rod_start_incline_angle=rod_start_incline_angle, + num_threads=num_threads, + save_data=True, + ) + + simulate_flow_past_rod_with_restart() diff --git a/sopht/utils/__init__.py b/sopht/utils/__init__.py index 8340ac4..d23d442 100644 --- a/sopht/utils/__init__.py +++ b/sopht/utils/__init__.py @@ -6,3 +6,4 @@ from .rod_viz import plot_video_of_rod_surface from .precision import get_real_t, get_test_tol from .pyst_kernel_config import get_pyst_dtype, get_pyst_kernel_config +from .restart_sim import restart_simulation diff --git a/sopht/utils/restart_sim.py b/sopht/utils/restart_sim.py new file mode 100644 index 0000000..f88de1d --- /dev/null +++ b/sopht/utils/restart_sim.py @@ -0,0 +1,44 @@ +from typing import Type +import os +import elastica as ea +from sopht.simulator.flow.navier_stokes_flow_simulators import ( + UnboundedNavierStokesFlowSimulator3D, + UnboundedNavierStokesFlowSimulator2D, +) +from sopht.simulator.flow.passive_transport_flow_simulators import ( + PassiveTransportFlowSimulator, +) +from sopht.utils.io import IO + + +def restart_simulation( + flow_sim: UnboundedNavierStokesFlowSimulator2D + | UnboundedNavierStokesFlowSimulator3D + | PassiveTransportFlowSimulator, + restart_simulator: Type[ea.BaseSystemCollection], + io: IO, + rod_io: IO, + forcing_io: IO, + restart_dir: str, +) -> None: + # find latest saved data + iter_num = [] + for filename in os.listdir(): + if "sopht" in filename and "h5" in filename: + iter_num.append(int(filename.split("_")[-1].split(".")[0])) + + if len(iter_num) == 0: + raise FileNotFoundError("There is no file to load in the directory.") + + latest = max(iter_num) + # load sopht data + flow_sim.time = io.load(h5_file_name=f"sopht_{latest:04d}.h5") + rod_io.load(h5_file_name=f"rod_{latest:04d}.h5") + forcing_io.load(h5_file_name=f"forcing_grid_{latest:04d}.h5") + rod_time = ea.load_state(restart_simulator, restart_dir, True) + + if flow_sim.time != rod_time: + raise ValueError( + "Simulation time of the flow is not matched with the Elastica, check your inputs!" + ) + print(f"sopht_{latest:04d}.h5 has been loaded") diff --git a/tests/test_utils/test_restart.py b/tests/test_utils/test_restart.py new file mode 100644 index 0000000..a2d0b62 --- /dev/null +++ b/tests/test_utils/test_restart.py @@ -0,0 +1,245 @@ +import pytest +import sopht.utils as spu +import sopht.simulator as sps +import elastica as ea +import numpy as np +import os +from numpy.testing import assert_allclose + + +@pytest.mark.parametrize("grid_size_x", [64]) +@pytest.mark.parametrize("precision", ["single"]) +@pytest.mark.parametrize("with_free_stream", [True]) +@pytest.mark.parametrize("filter_vorticity", [True]) +def test_restart_simulation(precision, grid_size_x, with_free_stream, filter_vorticity): + + final_time = 0.25 + + # run first half of the simulation + _ = run_sim( + final_time / 2, + save_data=True, + load_data=False, + precision=precision, + grid_size_x=grid_size_x, + with_free_stream=with_free_stream, + filter_vorticity=filter_vorticity, + ) + + # run second half of the simulation + recorded_data_restarted = run_sim( + final_time, + save_data=False, + load_data=True, + precision=precision, + grid_size_x=grid_size_x, + with_free_stream=with_free_stream, + filter_vorticity=filter_vorticity, + ) + + # run full simulation + recorded_data_full = run_sim( + final_time, + save_data=False, + load_data=False, + precision=precision, + grid_size_x=grid_size_x, + with_free_stream=with_free_stream, + filter_vorticity=filter_vorticity, + ) + + os.system("rm -f *h5 *xmf") + + for i in range(len(recorded_data_restarted)): + assert_allclose(recorded_data_restarted[i], recorded_data_full[i]) + + +def run_sim( + final_time: float, + save_data: bool, + load_data: bool, + precision: str, + grid_size_x: int, + with_free_stream: bool, + filter_vorticity: bool, +): + class RestartTestSimulator( + ea.BaseSystemCollection, ea.Constraints, ea.Forcing, ea.Damping + ): + ... + + restart_test_simulator = RestartTestSimulator() + + num_threads = 4 + grid_dim = 3 + x_range = 1.8 + nu = 3e-3 + real_t = spu.get_real_t(precision) + grid_size = (grid_size_x, int(0.25 * grid_size_x), grid_size_x) + filter_type = "multiplicative" + filter_order = 1 + + flow_sim = sps.UnboundedNavierStokesFlowSimulator3D( + grid_size=grid_size, + x_range=x_range, + kinematic_viscosity=nu, + with_forcing=True, + with_free_stream_flow=with_free_stream, + real_t=real_t, + num_threads=num_threads, + filter_vorticity=filter_vorticity, + filter_setting_dict={"type": filter_type, "order": filter_order}, + ) + + n_elems = 40 + cosserat_rod = ea.CosseratRod.straight_rod( + n_elems, + start=np.array([0.2 * x_range, 0.5 * x_range * 0.25, 0.75 * x_range]), + direction=np.array([0.0, 0.0, -1.0]), + normal=np.array([0.0, 1.0, 0.0]), + base_length=1, + base_radius=0.045, + density=830, + youngs_modulus=865, + shear_modulus=865 / 1.01, + ) + + restart_test_simulator.append(cosserat_rod) + + restart_test_simulator.constrain(cosserat_rod).using( + ea.OneEndFixedBC, + constrained_position_idx=(0,), + constrained_director_idx=(0,), + ) + # Add gravitational forces + restart_test_simulator.add_forcing_to(cosserat_rod).using( + ea.GravityForces, acc_gravity=np.array([0.0, 0.0, -0.036]) + ) + # add damping + dl = 1 / n_elems + rod_dt = 0.01 * dl + damping_constant = 1e-3 + restart_test_simulator.dampen(cosserat_rod).using( + ea.AnalyticalLinearDamper, + damping_constant=damping_constant, + time_step=rod_dt, + ) + + forcing_grid_cls = sps.CosseratRodSurfaceForcingGrid + rod_flow_interactor = sps.CosseratRodFlowInteraction( + cosserat_rod=cosserat_rod, + eul_grid_forcing_field=flow_sim.eul_grid_forcing_field, + eul_grid_velocity_field=flow_sim.velocity_field, + virtual_boundary_stiffness_coeff=-200000.0, + virtual_boundary_damping_coeff=-100.0, + dx=flow_sim.dx, + grid_dim=grid_dim, + real_t=real_t, + forcing_grid_cls=forcing_grid_cls, + surface_grid_density_for_largest_element=16, + ) + + restart_test_simulator.add_forcing_to(cosserat_rod).using( + sps.FlowForces, + rod_flow_interactor, + ) + + restart_test_simulator.finalize() + + restart_dir = "restart_data" + free_stream_velocity = np.array([1.0, 0.0, 0.0]) + real_t = spu.get_real_t(precision) + timestepper = ea.PositionVerlet() + do_step, stages_and_updates = ea.extend_stepper_interface( + timestepper, restart_test_simulator + ) + + recorded_data = [] + + # setup flow IO + io = spu.EulerianFieldIO( + position_field=flow_sim.position_field, + eulerian_fields_dict={ + "vorticity": flow_sim.vorticity_field, + "velocity": flow_sim.velocity_field, + }, + ) + # Initialize sphere IO + rod_io = spu.CosseratRodIO( + cosserat_rod=cosserat_rod, dim=grid_dim, real_dtype=real_t + ) + # Initialize forcing io + forcing_io = spu.IO(dim=grid_dim, real_dtype=real_t) + # Add vector field on lagrangian grid + forcing_io.add_as_lagrangian_fields_for_io( + lagrangian_grid=rod_flow_interactor.forcing_grid.position_field, + lagrangian_grid_name="cosseratrod", + vector_3d=rod_flow_interactor.lag_grid_forcing_field, + position_mismatch=rod_flow_interactor.lag_grid_position_mismatch_field, + ) + + if load_data: + spu.restart_simulation( + flow_sim=flow_sim, + restart_simulator=restart_test_simulator, + io=io, + rod_io=rod_io, + forcing_io=forcing_io, + restart_dir=restart_dir, + ) + + while flow_sim.time < final_time: + + if save_data: + io.save( + h5_file_name="sopht_" + str("%0.4d" % (flow_sim.time * 100)) + ".h5", + time=flow_sim.time, + ) + rod_io.save( + h5_file_name="rod_" + str("%0.4d" % (flow_sim.time * 100)) + ".h5", + time=flow_sim.time, + ) + forcing_io.save( + h5_file_name="forcing_grid_" + + str("%0.4d" % (flow_sim.time * 100)) + + ".h5", + time=flow_sim.time, + ) + ea.save_state(restart_test_simulator, restart_dir, flow_sim.time) + + # compute timestep + rod_dt = 2.5e-4 + flow_dt = rod_dt + + # timestep the rod, through the flow timestep + rod_time_steps = int(flow_dt / min(flow_dt, rod_dt)) + local_rod_dt = flow_dt / rod_time_steps + rod_time = flow_sim.time + for i in range(rod_time_steps): + rod_time = do_step( + timestepper, + stages_and_updates, + restart_test_simulator, + rod_time, + local_rod_dt, + ) + # timestep the cosserat_rod_flow_interactor + rod_flow_interactor.time_step(dt=local_rod_dt) + # evaluate feedback/interaction between flow and rod + rod_flow_interactor() + + flow_sim.time_step( + dt=flow_dt, + free_stream_velocity=free_stream_velocity, + ) + + recorded_data.append(flow_sim.velocity_field) + recorded_data.append(flow_sim.vorticity_field) + recorded_data.append(flow_sim.position_field) + recorded_data.append(np.array([flow_sim.time])) + recorded_data.append(rod_flow_interactor.lag_grid_forcing_field) + recorded_data.append(rod_flow_interactor.forcing_grid.position_field) + recorded_data.append(cosserat_rod.position_collection) + recorded_data.append(cosserat_rod.velocity_collection) + + return recorded_data