From 90adc4477c7ef469d0afe9190f0615daf3f7baa9 Mon Sep 17 00:00:00 2001 From: ilianasiriziba Date: Tue, 8 Aug 2023 11:50:23 -0500 Subject: [PATCH 01/11] restart function added --- sopht/utils/__init__.py | 1 + 1 file changed, 1 insertion(+) 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 From 91d43a914f7810b8edc9b86097c9bfdf8c6ce097 Mon Sep 17 00:00:00 2001 From: ilianasiriziba Date: Tue, 8 Aug 2023 11:51:17 -0500 Subject: [PATCH 02/11] initial commit of test for restart function --- tests/test_utils/test_restart.py | 255 +++++++++++++++++++++++++++++++ 1 file changed, 255 insertions(+) create mode 100644 tests/test_utils/test_restart.py diff --git a/tests/test_utils/test_restart.py b/tests/test_utils/test_restart.py new file mode 100644 index 0000000..f07cd02 --- /dev/null +++ b/tests/test_utils/test_restart.py @@ -0,0 +1,255 @@ +from typing import Type +import pytest +import sopht.utils as spu +import sopht.simulator as sps +import elastica as ea +import numpy as np +import os +from sopht.simulator.flow.flow_simulators import FlowSimulator +from sopht.simulator.immersed_body.immersed_body_flow_interaction import ( + ImmersedBodyFlowInteraction, +) +from numpy.testing import assert_allclose + + +@pytest.mark.parametrize("grid_size_x", [32]) +@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): + class RestartTestSimulator( + ea.BaseSystemCollection, ea.Constraints, ea.Forcing, ea.Damping + ): + ... + + for i in range(3): + restart_test_simulator = RestartTestSimulator() + + final_time = 0.25 + 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() + + if i == 0: + recorded_data = run_sim( + final_time / 2, + flow_sim, + restart_test_simulator, + rod_flow_interactor, + cosserat_rod, + save_data=True, + load_data=False, + grid_dim=grid_dim, + precision=precision, + ) + elif i == 1: + recorded_data = run_sim( + final_time, + flow_sim, + restart_test_simulator, + rod_flow_interactor, + cosserat_rod, + save_data=False, + load_data=True, + grid_dim=grid_dim, + precision=precision, + ) + else: + recorded_data_full = run_sim( + final_time, + flow_sim, + restart_test_simulator, + rod_flow_interactor, + cosserat_rod, + save_data=False, + load_data=False, + grid_dim=grid_dim, + precision=precision, + ) + + os.system("rm -f *h5 *xmf") + + for i in range(len(recorded_data)): + assert_allclose(recorded_data[i], recorded_data_full[i]) + + +def run_sim( + final_time: float, + flow_sim: Type[FlowSimulator], + restart_test_simulator: Type[ea.BaseSystemCollection], + rod_flow_interactor: Type[ImmersedBodyFlowInteraction], + cosserat_rod: Type[ea.CosseratRod], + save_data: bool, + load_data: bool, + grid_dim: int, + precision: str, +): + + 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=f"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(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 From 63b55e8aa46a0444e2beb83dc7c3411fbd9521b7 Mon Sep 17 00:00:00 2001 From: ilianasiriziba Date: Tue, 8 Aug 2023 11:51:49 -0500 Subject: [PATCH 03/11] initial commit of restart function --- sopht/utils/restart_sim.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 sopht/utils/restart_sim.py diff --git a/sopht/utils/restart_sim.py b/sopht/utils/restart_sim.py new file mode 100644 index 0000000..cc51ce7 --- /dev/null +++ b/sopht/utils/restart_sim.py @@ -0,0 +1,36 @@ +from typing import Type +import os +import elastica as ea +from sopht.simulator.flow.flow_simulators import FlowSimulator +from sopht.utils.io import IO + + +def restart_simulation( + flow_sim: Type[FlowSimulator], + 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: + print("There is no file to load in the directory") + return + + 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) + + assert ( + flow_sim.time == rod_time + ), "Simulation time of the flow is not matched with the Elastica, check your inputs!" + print(f"sopht_{latest:04d}.h5 has been loaded") From ac3fbacbe63ee4e7f617af7564892b4f0558d768 Mon Sep 17 00:00:00 2001 From: ilianasiriziba Date: Tue, 8 Aug 2023 11:52:09 -0500 Subject: [PATCH 04/11] initial commit of restart example --- .../RestartExample/restart_example.py | 425 ++++++++++++++++++ 1 file changed, 425 insertions(+) create mode 100644 examples/3d_examples/RestartExample/restart_example.py diff --git a/examples/3d_examples/RestartExample/restart_example.py b/examples/3d_examples/RestartExample/restart_example.py new file mode 100644 index 0000000..d0a1867 --- /dev/null +++ b/examples/3d_examples/RestartExample/restart_example.py @@ -0,0 +1,425 @@ +import elastica as ea +import numpy as np +import sopht.simulator as sps +import sopht.utils as spu +import click +import os +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", + vector_3d=cosserat_rod_flow_interactor.lag_grid_forcing_field, + position_mismatch=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() From b17e8c33773c627e8b157914ca2d7c1c51aa5a62 Mon Sep 17 00:00:00 2001 From: ilianasiriziba Date: Tue, 8 Aug 2023 12:03:58 -0500 Subject: [PATCH 05/11] initial commit of test for restart function --- tests/test_utils/test_restart.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_utils/test_restart.py b/tests/test_utils/test_restart.py index f07cd02..ea0213a 100644 --- a/tests/test_utils/test_restart.py +++ b/tests/test_utils/test_restart.py @@ -210,7 +210,7 @@ def run_sim( time=flow_sim.time, ) forcing_io.save( - h5_file_name=f"forcing_grid_" + h5_file_name="forcing_grid_" + str("%0.4d" % (flow_sim.time * 100)) + ".h5", time=flow_sim.time, From d5710b62baff87d8f019994f899c2689e231d74b Mon Sep 17 00:00:00 2001 From: ilianasiriziba Date: Tue, 8 Aug 2023 14:37:42 -0500 Subject: [PATCH 06/11] initial commit of restart function --- sopht/utils/restart_sim.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/sopht/utils/restart_sim.py b/sopht/utils/restart_sim.py index cc51ce7..9e6d776 100644 --- a/sopht/utils/restart_sim.py +++ b/sopht/utils/restart_sim.py @@ -1,12 +1,22 @@ -from typing import Type +from typing import Type, Union import os import elastica as ea -from sopht.simulator.flow.flow_simulators import FlowSimulator +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: Type[FlowSimulator], + flow_sim: Union[ + UnboundedNavierStokesFlowSimulator2D, + UnboundedNavierStokesFlowSimulator3D, + PassiveTransportFlowSimulator, + ], restart_simulator: Type[ea.BaseSystemCollection], io: IO, rod_io: IO, From e954142e5c5bdd2fe4f4bfb16feec6886f4ff6d5 Mon Sep 17 00:00:00 2001 From: ilianasiriziba Date: Tue, 8 Aug 2023 14:38:25 -0500 Subject: [PATCH 07/11] initial commit of test for restart function --- tests/test_utils/test_restart.py | 252 +++++++++++++++---------------- 1 file changed, 121 insertions(+), 131 deletions(-) diff --git a/tests/test_utils/test_restart.py b/tests/test_utils/test_restart.py index ea0213a..a2d0b62 100644 --- a/tests/test_utils/test_restart.py +++ b/tests/test_utils/test_restart.py @@ -1,160 +1,150 @@ -from typing import Type import pytest import sopht.utils as spu import sopht.simulator as sps import elastica as ea import numpy as np import os -from sopht.simulator.flow.flow_simulators import FlowSimulator -from sopht.simulator.immersed_body.immersed_body_flow_interaction import ( - ImmersedBodyFlowInteraction, -) from numpy.testing import assert_allclose -@pytest.mark.parametrize("grid_size_x", [32]) +@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): - class RestartTestSimulator( - ea.BaseSystemCollection, ea.Constraints, ea.Forcing, ea.Damping - ): - ... - - for i in range(3): - restart_test_simulator = RestartTestSimulator() - - final_time = 0.25 - 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, - ) + 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, + ) - restart_test_simulator.finalize() + # 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, + ) - if i == 0: - recorded_data = run_sim( - final_time / 2, - flow_sim, - restart_test_simulator, - rod_flow_interactor, - cosserat_rod, - save_data=True, - load_data=False, - grid_dim=grid_dim, - precision=precision, - ) - elif i == 1: - recorded_data = run_sim( - final_time, - flow_sim, - restart_test_simulator, - rod_flow_interactor, - cosserat_rod, - save_data=False, - load_data=True, - grid_dim=grid_dim, - precision=precision, - ) - else: - recorded_data_full = run_sim( - final_time, - flow_sim, - restart_test_simulator, - rod_flow_interactor, - cosserat_rod, - save_data=False, - load_data=False, - grid_dim=grid_dim, - precision=precision, - ) + # 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)): - assert_allclose(recorded_data[i], recorded_data_full[i]) + for i in range(len(recorded_data_restarted)): + assert_allclose(recorded_data_restarted[i], recorded_data_full[i]) def run_sim( final_time: float, - flow_sim: Type[FlowSimulator], - restart_test_simulator: Type[ea.BaseSystemCollection], - rod_flow_interactor: Type[ImmersedBodyFlowInteraction], - cosserat_rod: Type[ea.CosseratRod], save_data: bool, load_data: bool, - grid_dim: int, 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]) @@ -246,7 +236,7 @@ def run_sim( recorded_data.append(flow_sim.velocity_field) recorded_data.append(flow_sim.vorticity_field) recorded_data.append(flow_sim.position_field) - recorded_data.append(flow_sim.time) + 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) From 76d70c50f0f979602d8bc403fb6210d45b9a83d1 Mon Sep 17 00:00:00 2001 From: ilianasiriziba Date: Wed, 9 Aug 2023 11:58:58 -0500 Subject: [PATCH 08/11] initial commit of restart example --- examples/3d_examples/RestartExample/restart_example.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/3d_examples/RestartExample/restart_example.py b/examples/3d_examples/RestartExample/restart_example.py index d0a1867..b628d9d 100644 --- a/examples/3d_examples/RestartExample/restart_example.py +++ b/examples/3d_examples/RestartExample/restart_example.py @@ -157,8 +157,8 @@ class RestartExampleSimulator( forcing_io.add_as_lagrangian_fields_for_io( lagrangian_grid=cosserat_rod_flow_interactor.forcing_grid.position_field, lagrangian_grid_name="cosseratrod", - vector_3d=cosserat_rod_flow_interactor.lag_grid_forcing_field, - position_mismatch=cosserat_rod_flow_interactor.lag_grid_position_mismatch_field, + 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" From b4ee85dca5401c649043c543ea12924c47aa0c11 Mon Sep 17 00:00:00 2001 From: ilianasiriziba Date: Thu, 31 Aug 2023 16:34:06 -0500 Subject: [PATCH 09/11] Update: Union removed --- examples/3d_examples/RestartExample/restart_example.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/3d_examples/RestartExample/restart_example.py b/examples/3d_examples/RestartExample/restart_example.py index b628d9d..0189ea2 100644 --- a/examples/3d_examples/RestartExample/restart_example.py +++ b/examples/3d_examples/RestartExample/restart_example.py @@ -1,9 +1,8 @@ import elastica as ea import numpy as np -import sopht.simulator as sps import sopht.utils as spu +import sopht.simulator as sps import click -import os from matplotlib import pyplot as plt From f9faf00ddb73c4f8f739df9401fc8c6d7a245f32 Mon Sep 17 00:00:00 2001 From: ilianasiriziba Date: Thu, 31 Aug 2023 16:34:40 -0500 Subject: [PATCH 10/11] Update: Union removed --- sopht/utils/restart_sim.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/sopht/utils/restart_sim.py b/sopht/utils/restart_sim.py index 9e6d776..30d53a0 100644 --- a/sopht/utils/restart_sim.py +++ b/sopht/utils/restart_sim.py @@ -1,4 +1,4 @@ -from typing import Type, Union +from typing import Type import os import elastica as ea from sopht.simulator.flow.navier_stokes_flow_simulators import ( @@ -12,11 +12,9 @@ def restart_simulation( - flow_sim: Union[ - UnboundedNavierStokesFlowSimulator2D, - UnboundedNavierStokesFlowSimulator3D, - PassiveTransportFlowSimulator, - ], + flow_sim: UnboundedNavierStokesFlowSimulator2D + | UnboundedNavierStokesFlowSimulator3D + | PassiveTransportFlowSimulator, restart_simulator: Type[ea.BaseSystemCollection], io: IO, rod_io: IO, @@ -36,8 +34,8 @@ def restart_simulation( 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_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) assert ( From 962c756c1fb5cf74c14cfb4e6138e0066fc19009 Mon Sep 17 00:00:00 2001 From: ilianasiriziba Date: Tue, 5 Sep 2023 10:49:00 -0500 Subject: [PATCH 11/11] update: error raises if file not loaded correctly --- sopht/utils/restart_sim.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/sopht/utils/restart_sim.py b/sopht/utils/restart_sim.py index 30d53a0..f88de1d 100644 --- a/sopht/utils/restart_sim.py +++ b/sopht/utils/restart_sim.py @@ -28,8 +28,7 @@ def restart_simulation( iter_num.append(int(filename.split("_")[-1].split(".")[0])) if len(iter_num) == 0: - print("There is no file to load in the directory") - return + raise FileNotFoundError("There is no file to load in the directory.") latest = max(iter_num) # load sopht data @@ -38,7 +37,8 @@ def restart_simulation( forcing_io.load(h5_file_name=f"forcing_grid_{latest:04d}.h5") rod_time = ea.load_state(restart_simulator, restart_dir, True) - assert ( - flow_sim.time == rod_time - ), "Simulation time of the flow is not matched with the Elastica, check your inputs!" + 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")