|
| 1 | +#!/usr/bin/env python |
| 2 | +# encoding: utf-8 |
| 3 | +r""" |
| 4 | +Shock-tube problem |
| 5 | +=================================== |
| 6 | +
|
| 7 | +Solve the one-dimensional Euler equations for inviscid, compressible flow: |
| 8 | +
|
| 9 | +.. math:: |
| 10 | + \rho_t + (\rho u)_x & = 0 \\ |
| 11 | + (\rho u)_t + (\rho u^2 + p)_x & = 0 \\ |
| 12 | + E_t + (u (E + p) )_x & = 0. |
| 13 | +
|
| 14 | +The fluid is an ideal gas, with pressure given by :math:`p=\rho (\gamma-1)e` where |
| 15 | +e is internal energy. |
| 16 | +
|
| 17 | +This script runs a shock-tube problem. |
| 18 | +""" |
| 19 | +from clawpack import riemann |
| 20 | +from clawpack.riemann.euler_with_efix_1D_constants import * |
| 21 | + |
| 22 | +gamma = 1.4 # Ratio of specific heats |
| 23 | + |
| 24 | +def setup(use_petsc=False, outdir='./_output', solver_type='classic', |
| 25 | + kernel_language='Python',disable_output=False): |
| 26 | + |
| 27 | + if use_petsc: |
| 28 | + import clawpack.petclaw as pyclaw |
| 29 | + else: |
| 30 | + from clawpack import pyclaw |
| 31 | + |
| 32 | + if kernel_language =='Python': |
| 33 | + rs = riemann.euler_1D_py.euler_hllc_1D |
| 34 | + elif kernel_language =='Fortran': |
| 35 | + rs = riemann.euler_with_efix_1D |
| 36 | + |
| 37 | + if solver_type=='sharpclaw': |
| 38 | + solver = pyclaw.SharpClawSolver1D(rs) |
| 39 | + elif solver_type=='classic': |
| 40 | + solver = pyclaw.ClawSolver1D(rs) |
| 41 | + |
| 42 | + solver.kernel_language = kernel_language |
| 43 | + |
| 44 | + solver.bc_lower[0]=pyclaw.BC.extrap |
| 45 | + solver.bc_upper[0]=pyclaw.BC.extrap |
| 46 | + |
| 47 | + mx = 800; |
| 48 | + x = pyclaw.Dimension(-1.0,1.0,mx,name='x') |
| 49 | + domain = pyclaw.Domain([x]) |
| 50 | + state = pyclaw.State(domain,num_eqn) |
| 51 | + |
| 52 | + state.problem_data['gamma'] = gamma |
| 53 | + state.problem_data['gamma1'] = gamma - 1. |
| 54 | + |
| 55 | + x = state.grid.x.centers |
| 56 | + |
| 57 | + rho_l = 1.; rho_r = 1./8 |
| 58 | + p_l = 1.; p_r = 0.1 |
| 59 | + state.q[density ,:] = (x<0.)*rho_l + (x>=0.)*rho_r |
| 60 | + state.q[momentum,:] = 0. |
| 61 | + velocity = state.q[momentum,:]/state.q[density,:] |
| 62 | + pressure = (x<0.)*p_l + (x>=0.)*p_r |
| 63 | + state.q[energy ,:] = pressure/(gamma - 1.) + 0.5 * state.q[density,:] * velocity**2 |
| 64 | + |
| 65 | + claw = pyclaw.Controller() |
| 66 | + claw.tfinal = 0.4 |
| 67 | + claw.solution = pyclaw.Solution(state,domain) |
| 68 | + claw.solver = solver |
| 69 | + claw.num_output_times = 10 |
| 70 | + claw.outdir = outdir |
| 71 | + claw.setplot = setplot |
| 72 | + claw.keep_copy = True |
| 73 | + if disable_output: |
| 74 | + claw.output_format = None |
| 75 | + |
| 76 | + return claw |
| 77 | + |
| 78 | +#-------------------------- |
| 79 | +def setplot(plotdata): |
| 80 | +#-------------------------- |
| 81 | + """ |
| 82 | + Specify what is to be plotted at each frame. |
| 83 | + Input: plotdata, an instance of visclaw.data.ClawPlotData. |
| 84 | + Output: a modified version of plotdata. |
| 85 | + """ |
| 86 | + plotdata.clearfigures() # clear any old figures,axes,items data |
| 87 | + |
| 88 | + plotfigure = plotdata.new_plotfigure(name='', figno=0) |
| 89 | + |
| 90 | + plotaxes = plotfigure.new_plotaxes() |
| 91 | + plotaxes.axescmd = 'subplot(211)' |
| 92 | + plotaxes.title = 'Density' |
| 93 | + |
| 94 | + plotitem = plotaxes.new_plotitem(plot_type='1d') |
| 95 | + plotitem.plot_var = density |
| 96 | + plotitem.kwargs = {'linewidth':3} |
| 97 | + |
| 98 | + plotaxes = plotfigure.new_plotaxes() |
| 99 | + plotaxes.axescmd = 'subplot(212)' |
| 100 | + plotaxes.title = 'Energy' |
| 101 | + |
| 102 | + plotitem = plotaxes.new_plotitem(plot_type='1d') |
| 103 | + plotitem.plot_var = energy |
| 104 | + plotitem.kwargs = {'linewidth':3} |
| 105 | + |
| 106 | + return plotdata |
| 107 | + |
| 108 | +if __name__=="__main__": |
| 109 | + from clawpack.pyclaw.util import run_app_from_main |
| 110 | + output = run_app_from_main(setup,setplot) |
0 commit comments