1D Sod tube is a classic way to validate solvers for gas dynamic problems.

It looks like the following image. There is a 3D tube, which is shown as the blue wireframe, a diaphragm, which is shown as the red wire frame. The diaphragm divides the tube into two zones, which have different physical status of ideal gas separately. Please try to use show_diaphragm checkbox to show or hide the diaphragm. The bound slider could help you to zoom in and out. The tube is endless.

Now we are going to have a 1D tube. Let's try to make the tube radius apporach to zero, and then the tube is getting to be an 1D tube. Move the tube_radius slider to get that.

To begin with, we make a mesh on the tube. The grid points of the mesh are the location the CESE method uses to input and output the status values. Please try to check show_mesh to show the spatial part of the mesh.



In [1]:

    
# the magic inline command is used for ipython notebook figure output
# to embed the image on the page directly.
%matplotlib inline



In [1]:

    
import sodtube1d.handler_plot as hp
hp.interact_with_mesh_physical_model()

Actually, CESE method uses a mesh with grid points both in time and space. The mesh could be shown as the following:



In [2]:

    
import sodtube1d.generator_mesh as gmesh
generator_mesh = gmesh.Mesher()
generator_mesh.gen_mesh()
generator_mesh.show_mesh_in_space_time_ipynb(bound_x=0.0250, bound_t=0.02)

Please note the horizontal axis shows the spatial local and the vertical axis shows the time.

Now let us focus on the grid points at t = 0. This means only the grids of the mesh in the space are shown. They are the same thing as the mesh shown in the tube image above.



In [3]:

    
generator_mesh.show_mesh_in_space_time_ipynb(highlight=True, bound_x=0.0250, bound_t=0.02)

Secondly, we initialize the gas status on the grid points for this issue. This means we specify the values to the status of the ideal gas at t = 0.



In [3]:

    
import sodtube1d.solver_cese as cese
solver_cese = cese.Solver()
solver_cese.init_gas_status()

Let's have a look of the initial gas status.



In [4]:

    
import sodtube1d.helper_plot as helper_plot
solver_cese.data.refresh_solution()
solution = list(solver_cese.data.solution)
helper_plot.show_gas_status(solution)

From the figure above, it shows two points, one at the lefthand side of the diagphram and the other at the right, are specified initialized value and ready to iterate. Let's begin to iterate.



In [5]:

    
solver_cese._data.it_pt_nb = 2
solver_cese._data.it_nb = 0
solver_cese.cal_cese_status_before_half_dt()
solver_cese.cal_cese_status_after_half_dt()
solver_cese.push_status_along_t()

Let's have a look of the iterated result after 1 iteration from the image below. In the image, time pass by half of delta t which could be defined by users.



In [6]:

    
solver_cese.data.refresh_solution()
solution = list(solver_cese.data.solution)
helper_plot.show_gas_status(solution)

Everything looks good. Let's have one more iteration, that is, step in one more half of delta t and take a look the gas status.



In [7]:

    
solver_cese._data.it_nb = 1
solver_cese.cal_cese_status_before_half_dt()
solver_cese.cal_cese_status_after_half_dt()
solver_cese.push_status_along_t()
solver_cese.data.refresh_solution()
solution = list(solver_cese.data.solution)
helper_plot.show_gas_status(solution)



In [8]:

    
solver_cese._data.it_nb = 2
solver_cese.cal_cese_status_before_half_dt()
solver_cese.cal_cese_status_after_half_dt()
solver_cese.push_status_along_t()
solver_cese.data.refresh_solution()
solution = list(solver_cese.data.solution)
helper_plot.show_gas_status(solution)



In [3]:

    
import sodtube1d.generator_mesh as gmesh
import sodtube1d.solver_analytic as analytic
import sodtube1d.solver_cese as cese
import sodtube1d.handler_data as hd
import sodtube1d.handler_plot as hp


analytic_solver = analytic.Solver()
solution_analytic = analytic_solver.get_analytic_solution(mesh)

solver_cese = cese.Solver()
solution_cese = solver_cese.get_cese_solution()

dm = hd.DataManager()
pm = hp.PlotManager()

pm.get_plot_solutions_fig_rho(solution_analytic, solution_cese, "analytic rho", "cese rho")
pm.get_plot_solutions_fig_v(solution_analytic, solution_cese, "analytic v", "cese v")
pm.get_plot_solutions_fig_p(solution_analytic, solution_cese, "analytic p", "cese p")
pm.show_solution_comparison()

Then the two mesh points derive the status of a mesh point after half dt



In [ ]:

The whole process will look like the following. It shows the mesh points are used when iterating



In [ ]:

Now we get the status after time t. Further more, the time evolution of the status could be illistrated as the following



In [ ]: