1D Lagrangian shock-capturing hydro code

Riemann solver (this is done)

(* define function giving Mach number as a function of z ≡ ((Δv)/c) ≥ 0 *)
Clear[mach];
mach[z_,γ_]:=((γ+1)/4)z+(Sqrt[(((((γ+1)/4)z))^2)+1])

(* define function giving ratio ((Subscript[c, c])/c) ≡ ((sound speed at contact interface)/(sound speed outside boundary)) as a function of z ≡ ((Δv)/c), for both -ve and +ve z *)
Clear[ccc];
ccc[z_,γ_]:=If[z≤0,1+((γ-1)/2)z,(Sqrt[1+((γ-1)/2)z(-z+2mach[z,γ])])]

(* define function giving ratio ((Subscript[p, c])/p) ≡ ((pressure at contact interface)/(pressure outside boundary)) given z ≡ ((Δv)/c), for both -ve and +ve z *)
Clear[pcp];
pcp[z_,γ_]:=If[z≤0,((ccc[z,γ])^((2γ)/(γ-1))),1+γ z mach[z,γ]]

(* ((Subscript[p, c])/p) ≥ 0 for z ≥ zmin *)
Clear[zmin];
zmin[γ_]:=-(2/(γ-1))

(* define function that plots ((Subscript[p, c])/p) as a function of z ≡ ((Δv)/c) *)
Clear[plotpcp];
plotpcp[γ_]:=Plot[((pcp[z,γ])^((γ-1)/(2γ))),{z,zmin[γ],8}]

(* plot ((Subscript[p, c])/p) as a function of z ≡ ((Δv)/c) *)
plotpcp[(5/3)]

(* solve for velocity (Subscript[v, c]) and pressure (Subscript[p, c]) at contact interface by requiring that (Subscript[p, c]) be same for left and right sides *)
Clear[vcpc];
vcpc[{vl_,cl_,pl_},{vr_,cr_,pr_},γ_]:=(
Clear[vc];
vc=vc/.FindRoot[pl pcp[((vl-vc)/cl),γ]==pr pcp[((vc-vr)/cr),γ],{vc,vr+cr zmin[γ],vl-cl zmin[γ]}];
{vc,pl pcp[((vl-vc)/cl),γ]}
)

Integrator

(* integrate equations 1 step *)
Clear[step];
SetAttributes[step,HoldAll];
step[position_,zonevelocity_,zoneenergy_,velocity_,pressure_,zonemass_,Δt_]:=(
position+=;
zonevelocity+=;
zoneenergy+=;
)

Define global parameters

(* define adiabatic index *)
Clear[γ];
γ=

(* define mass of each Lagrangian zone *)
Clear[zonemass];
zonemass=

(* factor in Courant condition for timestep *)
Clear[courant];
courant=

(* define number of zones *)
Clear[Nzone];
Nzone=

(* size of box *)
Clear[boxsize];
boxsize=Nzone

Set initial conditions (this is done)

(* initialize density to top-hat function *)
Clear[zonedensity];
zonedensity=Table[If[(Nzone/4)≤i≤((3Nzone)/4),10,1],{i,(1/2),Nzone-(1/2)}];

(* deduce postions at edges of zones from zone densities *)
Clear[initialposition];
Do[If[i==0,initialposition={0},AppendTo[initialposition,initialposition[[i]]+(zonemass/(zonedensity[[i]]))]],{i,0,Nzone}]

(* adjust positions to desired boxsize *)
initialposition=Take[initialposition (boxsize/(initialposition[[-1]]-initialposition[[1]])),{1,Nzone}];

(* initialize velocities at centers of zones to zero *)
Clear[initialzonevelocity];
initialzonevelocity=Table[0,{i,1,Nzone}];

(* initialize sound speed c ≡ (Sqrt[γ(p/ρ)])to constant *)
Clear[zonesoundspeed];
zonesoundspeed=Table[1,{i,1,Nzone}];

(* deduce initial (energy/mass) ≡ (1/2)(v^2) + (1/(γ-1))(p/ρ) *)
Clear[initialzoneenergy];
initialzoneenergy=(1/2)(initialzonevelocity^2)+(1/(γ(γ-1)))(zonesoundspeed^2);

Define functions giving derived properties of zones

(* zone position x is the mean of the positions of its left and right edges *)
Clear[zonepositionfunc];
zonepositionfunc[position_,boxsize_]:=

(* zone density ρ *)
Clear[zonedensityfunc];
zonedensityfunc[position_,zonemass_,boxsize_]:=

(* zone pressure p *)
Clear[zonepressurefunc];
zonepressurefunc[zonevelocity_,zoneenergy_,zonedensity_,γ_]:=

(* zone sound speed c ≡ (Sqrt[γ(p/ρ)])*)
Clear[zonesoundspeedfunc];
zonesoundspeedfunc[zonedensity_,zonepressure_,γ_]:=

(* zone (entropy/mass) s ≡ Log[((c^(2/(1-γ)))/ρ)] *)
Clear[zoneentropyfunc];
zoneentropyfunc[zonedensity_,zonesoundspeed_,γ_]:=Log[((zonesoundspeed^(2/(γ-1)))/zonedensity)]

Time step

(* time step must be less than the sound-crossing time across any zone *)
Clear[Δtfunc];
Δtfunc[position_,zonesoundspeed_,courant_,boxsize_]:=

Hydro code (this is done, except for choosing nstep; you are welcome to do your own)

(* define a function to update evolutionary record of a variable *)
Clear[update];
SetAttributes[update,HoldAll];
update[variables_,variable_]:=If[!ListQ[variables],variables={variable},AppendTo[variables,variable]]

Clear[Δt,nstep];
(* number of steps to take *)
nstep=;
(* clear variables *)
Clear[position,zoneposition,zonevelocity,zoneenergy,zonedensity,zonepressure,zonesoundspeed,velocity,pressure,Δt];
(* clear arrays containing complete evolutionary record of each variable *)
Clear[positions,zonepositions,zonevelocities,zoneenergies,zonedensities,zonepressures,zonesoundspeeds,velocities,pressures,Δts];
(* initialize *)
position=initialposition;
zonevelocity=initialzonevelocity;
zoneenergy=initialzoneenergy;
(* each step *)
For[istep=0,istep≤nstep,istep++,
zoneposition=zonepositionfunc[position,boxsize];
zonedensity=zonedensityfunc[position,zonemass,boxsize];
zonepressure=zonepressurefunc[zonevelocity,zoneenergy,zonedensity,γ];
zonesoundspeed=zonesoundspeedfunc[zonedensity,zonepressure,γ];
(* Riemann solver for velocity and pressure at edges of zones *)
velocity={};
pressure={};
Do[
im=Mod[i-1,Nzone,1];
{v,p}=vcpc[{zonevelocity[[im]],zonesoundspeed[[im]],zonepressure[[im]]},{zonevelocity[[i]],zonesoundspeed[[i]],zonepressure[[i]]},γ];
AppendTo[velocity,v];
AppendTo[pressure,p];
,{i,1,Nzone}];
(* update evolutionary record of each variable *)
update[positions,position];
update[zonepositions,zoneposition];
update[zonevelocities,zonevelocity];
update[zoneenergies,zoneenergy];
update[zonedensities,zonedensity];
update[zonepressures,zonepressure];
update[zonesoundspeeds,zonesoundspeed];
update[velocities,velocity];
update[pressures,pressure];
(* time step *)
Δt=Δtfunc[position,zonesoundspeed,courant,boxsize];
update[Δts,Δt];
(* integrate *)
If[istep<nstep,step[position,zonevelocity,zoneenergy,velocity,pressure,zonemass,Δt]];
]

Clear[zoneentropies];
zoneentropies=zoneentropyfunc[zonedensities,zonesoundspeeds,γ];

Plot (you are welcome to do your own)

(* define function to plot several things side by side *)
Clear[plots];
plots[plotlist_]:=GraphicsRow[plotlist,0,ImageSize->400Length[plotlist]]

(* define command to plot variables at centers of zones *)
Clear[plot]
plot[positions_,variables_,istep_,plotrange_]:=If[0≤istep≤nstep,ListLinePlot[Table[{positions[[istep+1,i]],variables[[istep+1,i]]},{i,1,Nzone}],PlotRange->plotrange,PlotMarkers->{Automatic,4},Frame->True,Axes->False]]

Animate[plot[zonepositions,zonepressures,istep,{0,3.4}],{istep,0,nstep,1},DefaultDuration->30,AnimationRepetitions->1,AnimationRunning->False]

following plots are for istep = 100, γ = (5/3), zonemass = 1, courant = 0.7, Nzone = (2^8), boxsize = Nzone

plots[{plot[zonepositions,zonepressures,100,{0,All}],plot[positions,pressures,100,{0,All}]}]

plots[{plot[zonepositions,zonevelocities,100,{All,All}],plot[positions,velocities,100,{All,All}]}]

plots[{plot[zonepositions,zonesoundspeeds,100,{0,All}],plot[zonepositions,zonedensities,100,{0,All}]}]