Merge pull request #104 from SophT-Team/103_cleanup_flow_simulator

(#103) cleanup flow simulator
SophT-Team · Jan 14, 2023 · a42e862 · a42e862
2 parents afe9861 + 712ad61
commit a42e862
Show file tree

Hide file tree

Showing 26 changed files with 1,548 additions and 1,097 deletions.
diff --git a/examples/2d_examples/FlowPastCylinderCase/flow_past_cylinder.py b/examples/2d_examples/FlowPastCylinderCase/flow_past_cylinder.py
@@ -32,18 +32,16 @@ def flow_past_cylinder_boundary_forcing_case(
     cyl_radius = 0.03
     nu = cyl_radius * velocity_scale / reynolds
     x_range = 1.0
-
-    flow_sim = sps.UnboundedFlowSimulator2D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator2D(
         grid_size=grid_size,
         x_range=x_range,
         kinematic_viscosity=nu,
-        flow_type="navier_stokes_with_forcing",
+        with_forcing=True,
         with_free_stream_flow=True,
         real_t=real_t,
         num_threads=num_threads,
         time=0.0,
     )
-
     # Initialize fixed cylinder (elastica rigid body) with direction along Z
     x_cm = 2.5 * cyl_radius
     y_cm = 0.5 * grid_size_y / grid_size_x

diff --git a/examples/2d_examples/FlowPastRodCase/flow_past_rod.py b/examples/2d_examples/FlowPastRodCase/flow_past_rod.py
@@ -101,11 +101,11 @@ class FlowPastRodSimulator(
     # Flow parameters
     # Re = velocity_free_stream * base_length / nu
     nu = base_length * velocity_free_stream / reynolds
-    flow_sim = sps.UnboundedFlowSimulator2D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator2D(
         grid_size=(grid_size_y, grid_size_x),
         x_range=x_range,
         kinematic_viscosity=nu,
-        flow_type="navier_stokes_with_forcing",
+        with_forcing=True,
         with_free_stream_flow=True,
         real_t=real_t,
         num_threads=num_threads,

diff --git a/examples/2d_examples/ImmersedContinuumSnakeCase/immersed_continuum_snake.py b/examples/2d_examples/ImmersedContinuumSnakeCase/immersed_continuum_snake.py
@@ -86,11 +86,11 @@ class ImmersedContinuumSnakeSimulator(
     # Flow parameters
     vel_scale = base_length / period
     nu = base_length * vel_scale / reynolds
-    flow_sim = sps.UnboundedFlowSimulator2D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator2D(
         grid_size=grid_size,
         x_range=x_range,
         kinematic_viscosity=nu,
-        flow_type="navier_stokes_with_forcing",
+        with_forcing=True,
         real_t=real_t,
         num_threads=num_threads,
     )

diff --git a/..._examples/ImmersedFlexiblePendulumCases/immersed_flexible_pendulum_with_rigid_cylinder.py b/..._examples/ImmersedFlexiblePendulumCases/immersed_flexible_pendulum_with_rigid_cylinder.py
@@ -91,11 +91,11 @@ class ImmersedFlexiblePendulumSimulator(
     vel_scale = np.sqrt(np.fabs(gravitational_acc) * rod_length)
     Re = 500
     nu = rod_length * vel_scale / Re
-    flow_sim = sps.UnboundedFlowSimulator2D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator2D(
         grid_size=grid_size,
         x_range=x_range,
         kinematic_viscosity=nu,
-        flow_type="navier_stokes_with_forcing",
+        with_forcing=True,
         real_t=real_t,
         num_threads=num_threads,
     )

diff --git a/examples/2d_examples/LambOseenVortexCase/lamb_oseen_vortex.py b/examples/2d_examples/LambOseenVortexCase/lamb_oseen_vortex.py
@@ -27,17 +27,16 @@ def lamb_oseen_vortex_flow_case(
     t_end = 1.4
     # to start with max circulation = 1
     gamma = 4 * np.pi * nu * t_start
-    flow_sim = sps.UnboundedFlowSimulator2D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator2D(
         grid_size=grid_size,
         x_range=x_range,
         kinematic_viscosity=nu,
-        flow_type="navier_stokes",
+        with_forcing=False,
         with_free_stream_flow=True,
         real_t=real_t,
         num_threads=num_threads,
         time=t_start,
     )
-
     flow_sim.vorticity_field[...] = compute_lamb_oseen_vorticity(
         x=flow_sim.position_field[x_axis_idx],
         y=flow_sim.position_field[y_axis_idx],
@@ -128,22 +127,6 @@ def lamb_oseen_vortex_flow_case(
     print(f"Final vortex center location: ({x_cm_final}, {y_cm_final})")
     print(f"vorticity L2 error: {l2_error}")
     print(f"vorticity Linf error: {linf_error}")
-    final_analytical_velocity_field = compute_lamb_oseen_velocity(
-        x=flow_sim.position_field[x_axis_idx],
-        y=flow_sim.position_field[y_axis_idx],
-        x_cm=x_cm_final,
-        y_cm=y_cm_final,
-        nu=nu,
-        gamma=gamma,
-        t=t_end,
-        real_t=real_t,
-    )
-    flow_sim.compute_velocity_from_vorticity()
-    error_field = np.fabs(flow_sim.velocity_field - final_analytical_velocity_field)
-    l2_error = np.linalg.norm(error_field) * flow_sim.dx
-    linf_error = np.amax(error_field)
-    print(f"velocity L2 error: {l2_error}")
-    print(f"velocity Linf error: {linf_error}")
 
 
 if __name__ == "__main__":

diff --git a/examples/2d_examples/LidDrivenCavityCase/lid_driven_cavity_case.py b/examples/2d_examples/LidDrivenCavityCase/lid_driven_cavity_case.py
@@ -27,15 +27,14 @@ def lid_driven_cavity_case(
     ldc_side_length = 0.6
     nu = ldc_side_length * lid_velocity / reynolds
     x_range = 1.0
-    flow_sim = sps.UnboundedFlowSimulator2D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator2D(
         grid_size=grid_size,
         x_range=x_range,
         kinematic_viscosity=nu,
-        flow_type="navier_stokes_with_forcing",
+        with_forcing=True,
         real_t=real_t,
         num_threads=num_threads,
     )
-
     # Initialize lid driven cavity forcing grid
     num_lag_nodes_per_side = grid_size[x_axis_idx] * 3 // 8
     ldc_com = (0.5, 0.5)

diff --git a/examples/2d_examples/OctopusArmCase/run_coupling_tapered_arm_and_cylinder.py b/examples/2d_examples/OctopusArmCase/run_coupling_tapered_arm_and_cylinder.py
@@ -66,11 +66,11 @@ def tapered_arm_and_cylinder_flow_coupling(
     # Flow parameters
     vel_scale = base_length / period
     nu = base_length * vel_scale / reynolds
-    flow_sim = sps.UnboundedFlowSimulator2D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator2D(
         grid_size=grid_size,
         x_range=x_range,
         kinematic_viscosity=nu,
-        flow_type="navier_stokes_with_forcing",
+        with_forcing=True,
         real_t=real_t,
         num_threads=num_threads,
     )

diff --git a/examples/2d_examples/VortexInducedFlagFlappingCase/turek_and_hron_benchmark.py b/examples/2d_examples/VortexInducedFlagFlappingCase/turek_and_hron_benchmark.py
@@ -99,11 +99,11 @@ class FlowPastRodSimulator(ea.BaseSystemCollection, ea.Constraints, ea.Forcing):
     # Re = velocity_free_stream * cyl_diameter / nu
     cyl_diameter = base_length * cyl_diameter_to_rod_length
     nu = cyl_diameter * velocity_free_stream / reynolds
-    flow_sim = sps.UnboundedFlowSimulator2D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator2D(
         grid_size=grid_size,
         x_range=x_range,
         kinematic_viscosity=nu,
-        flow_type="navier_stokes_with_forcing",
+        with_forcing=True,
         with_free_stream_flow=True,
         real_t=real_t,
         num_threads=num_threads,

diff --git a/examples/3d_examples/ElasticNetCase/immersed_elastic_net.py b/examples/3d_examples/ElasticNetCase/immersed_elastic_net.py
@@ -38,11 +38,11 @@ def immersed_elastic_net_case(
         * vel_free_stream
         / reynolds
     )
-    flow_sim = sps.UnboundedFlowSimulator3D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator3D(
         grid_size=grid_size,
         x_range=domain_x_range,
         kinematic_viscosity=kinematic_viscosity,
-        flow_type="navier_stokes_with_forcing",
+        with_forcing=True,
         with_free_stream_flow=True,
         real_t=real_t,
         num_threads=num_threads,

diff --git a/examples/3d_examples/FlowPastRodCase/flow_past_rod_case.py b/examples/3d_examples/FlowPastRodCase/flow_past_rod_case.py
@@ -99,11 +99,11 @@ class FlowPastRodSimulator(
 
     # ==================FLOW SETUP START=========================
     kinematic_viscosity = U_free_stream * base_diameter / reynolds
-    flow_sim = sps.UnboundedFlowSimulator3D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator3D(
         grid_size=grid_size,
         x_range=x_range,
         kinematic_viscosity=kinematic_viscosity,
-        flow_type="navier_stokes_with_forcing",
+        with_forcing=True,
         with_free_stream_flow=True,
         real_t=real_t,
         num_threads=num_threads,

diff --git a/examples/3d_examples/FlowPastSphereCase/flow_past_sphere_case.py b/examples/3d_examples/FlowPastSphereCase/flow_past_sphere_case.py
@@ -28,13 +28,13 @@ def flow_past_sphere_case(
     far_field_velocity = 1.0
     sphere_diameter = 0.4 * min(grid_size_z, grid_size_y) / grid_size_x * x_range
     nu = far_field_velocity * sphere_diameter / reynolds
-    flow_sim = sps.UnboundedFlowSimulator3D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator3D(
         grid_size=grid_size,
         x_range=x_range,
         kinematic_viscosity=nu,
         real_t=real_t,
         num_threads=num_threads,
-        flow_type="navier_stokes_with_forcing",
+        with_forcing=True,
         with_free_stream_flow=True,
         time=0.0,
         # filter_vorticity=True,

diff --git a/examples/3d_examples/FlowThroughPipeCase/flow_through_circular_pipe_case.py b/examples/3d_examples/FlowThroughPipeCase/flow_through_circular_pipe_case.py
@@ -29,11 +29,11 @@ def flow_through_circular_pipe_case(
     z_axis_idx = spu.VectorField.z_axis_idx()
     x_range = 1.0
     nu = 1e-2
-    flow_sim = sps.UnboundedFlowSimulator3D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator3D(
         grid_size=grid_size,
         x_range=x_range,
         kinematic_viscosity=nu,
-        flow_type="navier_stokes_with_forcing",
+        with_forcing=True,
         with_free_stream_flow=True,
         real_t=real_t,
         num_threads=num_threads,

diff --git a/examples/3d_examples/HillSphericalVortexCase/hill_sphere_vortex_case.py b/examples/3d_examples/HillSphericalVortexCase/hill_sphere_vortex_case.py
@@ -23,14 +23,12 @@ def hill_sphere_vortex_case(
     z_axis_idx = spu.VectorField.z_axis_idx()
     # Consider a 1 by 1 by 1 3D domain
     x_range = 1.0
-    flow_sim = sps.UnboundedFlowSimulator3D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator3D(
         grid_size=grid_size,
         x_range=x_range,
         kinematic_viscosity=0.0,
-        flow_type="navier_stokes",
         real_t=real_t,
         num_threads=num_threads,
-        navier_stokes_inertial_term_form="advection_stretching_split",
     )
     # init vortex at domain center
     vortex_origin = (
@@ -53,7 +51,17 @@ def hill_sphere_vortex_case(
         y_grid=flow_sim.position_field[y_axis_idx],
         z_grid=flow_sim.position_field[z_axis_idx],
     )
-    flow_sim.compute_flow_velocity(free_stream_velocity=np.zeros(grid_dim))
+    # compute velocity from vorticity
+    flow_sim._unbounded_poisson_solver.vector_field_solve(
+        solution_vector_field=flow_sim.stream_func_field,
+        rhs_vector_field=flow_sim.vorticity_field,
+    )
+    flow_sim._curl(
+        curl=flow_sim.velocity_field,
+        field=flow_sim.stream_func_field,
+        prefactor=flow_sim.real_t(0.5 / flow_sim.dx),
+    )
+    # compute diagnostics
     numerical_kinetic_energy = (
         0.5 * np.sum(np.square(flow_sim.velocity_field)) * flow_sim.dx**3
     )
@@ -107,13 +115,6 @@ def hill_sphere_vortex_case(
     ax2.set_ylabel("Radial velocity")
     fig2.savefig("midplane_radial_velocity.png")
 
-    # check the vorticity stretching term
-    _, _, _, _, sphere_r_grid = hill_sphere_vortex.compute_local_coordinates(
-        flow_sim.position_field[x_axis_idx],
-        flow_sim.position_field[y_axis_idx],
-        flow_sim.position_field[z_axis_idx],
-    )
-
     if save_data:
         # setup IO
         # TODO internalise this in flow simulator as dump_fields

diff --git a/examples/3d_examples/MagneticCiliaCarpetCase/immersed_magnetic_cilia_carpet.py b/examples/3d_examples/MagneticCiliaCarpetCase/immersed_magnetic_cilia_carpet.py
@@ -32,11 +32,11 @@ def immersed_magnetic_cilia_carpet_case(
         * cilia_carpet_simulator.velocity_scale
         / reynolds
     )
-    flow_sim = sps.UnboundedFlowSimulator3D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator3D(
         grid_size=grid_size,
         x_range=domain_x_range,
         kinematic_viscosity=kinematic_viscosity,
-        flow_type="navier_stokes_with_forcing",
+        with_forcing=True,
         real_t=real_t,
         num_threads=num_threads,
         filter_vorticity=True,

diff --git a/examples/3d_examples/OctopusArmCase3D/run_tapered_arm_and_sphere_with_flow.py b/examples/3d_examples/OctopusArmCase3D/run_tapered_arm_and_sphere_with_flow.py
@@ -126,11 +126,11 @@ def tapered_arm_and_cylinder_flow_coupling(
     # ==================FLOW SETUP START=========================
     # Flow parameters
     kinematic_viscosity = base_diameter * vel_scale / reynolds_number
-    flow_sim = sps.UnboundedFlowSimulator3D(
+    flow_sim = sps.UnboundedNavierStokesFlowSimulator3D(
         grid_size=grid_size,
         x_range=x_range,
         kinematic_viscosity=kinematic_viscosity,
-        flow_type="navier_stokes_with_forcing",
+        with_forcing=True,
         with_free_stream_flow=False,
         real_t=real_t,
         num_threads=num_threads,

diff --git a/examples/3d_examples/PointSourceAdvectAndDiffuseCase/point_source_advection_diffusion.py b/examples/3d_examples/PointSourceAdvectAndDiffuseCase/point_source_advection_diffusion.py
@@ -26,22 +26,23 @@ def point_source_advection_diffusion_case(
     # start with non-zero to avoid singularity in point source
     t_start = 5.0
     t_end = 5.4
-    # init vortex at (0.3 0.3, 0.3)
-    flow_sim = sps.UnboundedFlowSimulator3D(
+    flow_sim = sps.PassiveTransportFlowSimulator(
+        kinematic_viscosity=nu,
+        grid_dim=grid_dim,
         grid_size=grid_size,
         x_range=x_range,
-        kinematic_viscosity=nu,
-        flow_type="passive_vector",
         real_t=real_t,
         num_threads=num_threads,
         time=t_start,
+        field_type="vector",
     )
+    # init vortex at (0.3 0.3, 0.3)
     x_cm_start = 0.3
     y_cm_start = x_cm_start
     z_cm_start = x_cm_start
     # to start with point source magnitude = 1
     point_mag = 4.0 * np.pi * nu * t_start**1.5
-    vorticity_field = flow_sim.primary_vector_field.view()
+    vorticity_field = flow_sim.primary_field.view()
     for i in range(grid_dim):
         vorticity_field[i] = compute_diffused_point_source_field(
             x_grid=flow_sim.position_field[x_axis_idx],

diff --git a/sopht/simulator/flow/__init__.py b/sopht/simulator/flow/__init__.py
@@ -1,2 +1,8 @@
 from .flow_simulators_2d import UnboundedFlowSimulator2D
 from .flow_simulators_3d import UnboundedFlowSimulator3D
+from .flow_simulators import FlowSimulator
+from .passive_transport_flow_simulators import PassiveTransportFlowSimulator
+from .navier_stokes_flow_simulators import (
+    UnboundedNavierStokesFlowSimulator2D,
+    UnboundedNavierStokesFlowSimulator3D,
+)