Sod Shock Tube
The Sod shock tube is the classic test problem for compressible flow solvers. It consists of a 1D Riemann problem with a high-pressure region on the left and a low-pressure region on the right, separated by a diaphragm at $x = 0.5$. The resulting wave structure includes a left-moving rarefaction fan, a right-moving contact discontinuity, and a right-moving shock wave.
Problem Setup
The 1D Euler equations are:
\[\pdv{}{t}\begin{pmatrix}\rho \\ \rho v \\ E\end{pmatrix} + \pdv{}{x}\begin{pmatrix}\rho v \\ \rho v^2 + P \\ (E+P)v\end{pmatrix} = 0,\]
with initial conditions:
\[(\rho, v, P) = \begin{cases}(1, 0, 1) & x < 0.5,\\(0.125, 0, 0.1) & x \geq 0.5.\end{cases}\]
We begin by loading the package and defining the problem.
using FiniteVolumeMethod
using StaticArrays
gamma = 1.4
eos = IdealGasEOS(gamma)
law = EulerEquations{1}(eos)FiniteVolumeMethod.EulerEquations{1, FiniteVolumeMethod.IdealGasEOS{Float64}}(FiniteVolumeMethod.IdealGasEOS{Float64}(1.4))Define the left and right primitive states $(\rho, v, P)$:
wL = SVector(1.0, 0.0, 1.0)
wR = SVector(0.125, 0.0, 0.1)3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
0.125
0.0
0.1Set up the mesh, boundary conditions, and initial condition:
N = 200
mesh = StructuredMesh1D(0.0, 1.0, N)
bc_left = DirichletHyperbolicBC(wL)
bc_right = DirichletHyperbolicBC(wR)FiniteVolumeMethod.DirichletHyperbolicBC{3, Float64}([0.125, 0.0, 0.1])The initial condition places the diaphragm at $x = 0.5$:
ic(x) = x < 0.5 ? wL : wRic (generic function with 1 method)Solving with HLLC + MUSCL
We use the HLLC Riemann solver (which resolves the contact discontinuity) with MUSCL reconstruction using the minmod limiter for second-order accuracy.
prob_hllc = HyperbolicProblem(
law, mesh, HLLCSolver(), CellCenteredMUSCL(MinmodLimiter()),
bc_left, bc_right, ic;
final_time = 0.2, cfl = 0.5
)
x_hllc, U_hllc, t_hllc = solve_hyperbolic(prob_hllc)200-element Vector{Float64}:
0.0025
0.0075
0.0125
0.0175
0.0225
⋮
0.9775
0.9825
0.9875
0.9925
0.9975Solving with HLL + MUSCL
For comparison, we also solve with the HLL solver, which does not resolve the contact discontinuity and is therefore more diffusive.
prob_hll = HyperbolicProblem(
law, mesh, HLLSolver(), CellCenteredMUSCL(MinmodLimiter()),
bc_left, bc_right, ic;
final_time = 0.2, cfl = 0.5
)
x_hll, U_hll, t_hll = solve_hyperbolic(prob_hll)([0.0025, 0.0075, 0.0125, 0.0175, 0.0225, 0.0275, 0.0325, 0.0375, 0.0425, 0.0475, 0.0525, 0.0575, 0.0625, 0.0675, 0.0725, 0.0775, 0.0825, 0.08750000000000001, 0.0925, 0.0975, 0.10250000000000001, 0.1075, 0.1125, 0.11750000000000001, 0.1225, 0.1275, 0.1325, 0.1375, 0.14250000000000002, 0.1475, 0.1525, 0.1575, 0.1625, 0.1675, 0.17250000000000001, 0.1775, 0.1825, 0.1875, 0.1925, 0.1975, 0.2025, 0.20750000000000002, 0.2125, 0.2175, 0.2225, 0.2275, 0.2325, 0.23750000000000002, 0.2425, 0.2475, 0.2525, 0.2575, 0.2625, 0.2675, 0.2725, 0.2775, 0.28250000000000003, 0.28750000000000003, 0.2925, 0.2975, 0.3025, 0.3075, 0.3125, 0.3175, 0.3225, 0.3275, 0.3325, 0.3375, 0.3425, 0.34750000000000003, 0.3525, 0.3575, 0.3625, 0.3675, 0.3725, 0.3775, 0.3825, 0.3875, 0.3925, 0.3975, 0.4025, 0.40750000000000003, 0.41250000000000003, 0.4175, 0.4225, 0.4275, 0.4325, 0.4375, 0.4425, 0.4475, 0.4525, 0.4575, 0.4625, 0.4675, 0.47250000000000003, 0.47750000000000004, 0.4825, 0.4875, 0.4925, 0.4975, 0.5025000000000001, 0.5075000000000001, 0.5125, 0.5175, 0.5225, 0.5275, 0.5325, 0.5375, 0.5425, 0.5475, 0.5525, 0.5575, 0.5625, 0.5675, 0.5725, 0.5775, 0.5825, 0.5875, 0.5925, 0.5975, 0.6025, 0.6075, 0.6125, 0.6175, 0.6225, 0.6275000000000001, 0.6325000000000001, 0.6375000000000001, 0.6425, 0.6475, 0.6525, 0.6575, 0.6625, 0.6675, 0.6725, 0.6775, 0.6825, 0.6875, 0.6925, 0.6975, 0.7025, 0.7075, 0.7125, 0.7175, 0.7225, 0.7275, 0.7325, 0.7375, 0.7425, 0.7475, 0.7525000000000001, 0.7575000000000001, 0.7625000000000001, 0.7675000000000001, 0.7725, 0.7775, 0.7825, 0.7875, 0.7925, 0.7975, 0.8025, 0.8075, 0.8125, 0.8175, 0.8225, 0.8275, 0.8325, 0.8375, 0.8425, 0.8475, 0.8525, 0.8575, 0.8625, 0.8675, 0.8725, 0.8775000000000001, 0.8825000000000001, 0.8875000000000001, 0.8925000000000001, 0.8975, 0.9025, 0.9075, 0.9125, 0.9175, 0.9225, 0.9275, 0.9325, 0.9375, 0.9425, 0.9475, 0.9525, 0.9575, 0.9625, 0.9675, 0.9725, 0.9775, 0.9825, 0.9875, 0.9925, 0.9975], StaticArraysCore.SVector{3, Float64}[[1.0, 0.0, 2.5000000000000004], [1.0, 0.0, 2.5000000000000004], [1.0, 0.0, 2.5000000000000004], [1.0, 0.0, 2.5000000000000004], [1.0, 0.0, 2.5000000000000004], [1.0, 0.0, 2.5000000000000004], [1.0, 0.0, 2.5000000000000004], [1.0, 0.0, 2.5000000000000004], [1.0, 0.0, 2.5000000000000004], [1.0, 0.0, 2.5000000000000004], [1.0, 0.0, 2.5000000000000004], [1.0, 0.0, 2.5000000000000004], [1.0, 0.0, 2.5000000000000004], [1.0, 8.437252885200764e-18, 2.5000000000000004], [0.9999999999999994, 2.491279178948263e-16, 2.5], [0.9999999999999986, 1.213680291025257e-15, 2.499999999999997], [0.999999999999996, 4.1691767448704535e-15, 2.4999999999999885], [0.9999999999999883, 1.3078185760653913e-14, 2.4999999999999614], [0.9999999999999649, 4.09834271281666e-14, 2.499999999999879], [0.9999999999998919, 1.2720691572341375e-13, 2.4999999999996234], [0.9999999999996673, 3.9286134778513786e-13, 2.4999999999988365], [0.999999999998985, 1.1998131006198793e-12, 2.499999999996449], [0.9999999999969469, 3.6114628850721256e-12, 2.499999999989316], [0.9999999999909579, 1.069754186316371e-11, 2.499999999968355], [0.9999999999736584, 3.1166458400476345e-11, 2.499999999907807], [0.9999999999245529, 8.926852518971106e-11, 2.499999999735938], [0.999999999787645, 2.5125974546761897e-10, 2.499999999256761], [0.999999999412885, 6.946811740866965e-10, 2.4999999979451024], [0.9999999984061687, 1.885843744361201e-9, 2.4999999944215956], [0.9999999957535254, 5.024493835977588e-9, 2.499999985137343], [0.9999999889010945, 1.3132399460265825e-8, 2.499999961153835], [0.9999999715565502, 3.365474058657993e-8, 2.499999900447931], [0.9999999285659291, 8.452192501341613e-8, 2.4999997499807676], [0.9999998242871386, 2.0790623030559122e-7, 2.499999385005054], [0.9999995769300942, 5.005828990029401e-7, 2.499998519255697], [0.9999990035728858, 1.1789875812742161e-6, 2.4999965125070647], [0.9999977059953318, 2.7142983622277284e-6, 2.4999919709938143], [0.9999948415192289, 6.103574086344602e-6, 2.4999819453675], [0.9999886796310106, 1.3394333877536001e-5, 2.4999603789450964], [0.999975778112993, 2.8659239783625316e-5, 2.4999152244557767], [0.9999495197137271, 5.972701362669346e-5, 2.4998233235088536], [0.9998976433314871, 0.00012110167592520782, 2.499641769833816], [0.9997983221704386, 0.00023859641437283806, 2.499294196768471], [0.9996143727822068, 0.0004561648650034905, 2.4986505528076535], [0.9992854683038926, 0.0008450542651484602, 2.4974999749578184], [0.9987189515570964, 0.0015145154402055468, 2.4955189691460316], [0.9977810972343429, 0.002621758249274499, 2.4922416209489597], [0.9962922707553308, 0.004376864492871164, 2.487044318710204], [0.9940306702305028, 0.007036848970262733, 2.479161896073423], [0.9907489897104822, 0.010883568779716144, 2.467750644360798], [0.9862052304416464, 0.01618425862880229, 2.452001916791638], [0.9802030829550274, 0.023141022693239123, 2.4312883702215453], [0.9726311320569815, 0.03184351512560171, 2.4053022214778315], [0.9634875577119812, 0.04224195350990727, 2.374136379311961], [0.9528808606726294, 0.0541518179401192, 2.338275431220649], [0.9410062601586417, 0.06728904569607308, 2.2984987766517406], [0.9281067068555223, 0.08132249705397444, 2.2557324272555928], [0.9144316205063826, 0.09592601920878363, 2.210899245018311], [0.9002039323038724, 0.11081704846786278, 2.1648051559097876], [0.8856002892407482, 0.12577808777202643, 2.11807301491319], [0.8707547570835176, 0.1406644713266877, 2.0711210023476783], [0.8558145459582368, 0.15526014075403272, 2.024462889410166], [0.840919979392585, 0.16945003253041452, 1.97863168200543], [0.8261248585620014, 0.18318464273046264, 1.9337251428378468], [0.8114639788623139, 0.1964430649194128, 1.8898261242036143], [0.7969569385720822, 0.20920992720606937, 1.8469736285877414], [0.7826176900160137, 0.22147756209284525, 1.8051878881272658], [0.7684575634676222, 0.2332442239809894, 1.7644800273462091], [0.7544851480226102, 0.24451099669005152, 1.7248517702522705], [0.7407072048755996, 0.25528113224447435, 1.6862986569129588], [0.7271293669598613, 0.26555990939838914, 1.6488125248285372], [0.7137555198165886, 0.2753580325608214, 1.6123813695007625], [0.7005817328367857, 0.2846811826081359, 1.5769679125346725], [0.6876091362158797, 0.29354357006562304, 1.5425537652037868], [0.6748365774739121, 0.30195635404776333, 1.5091131074508382], [0.6622602802583001, 0.3099268003324023, 1.4766113070551163], [0.6498792041205113, 0.3174663808743362, 1.4450240867570652], [0.6376960473895368, 0.3245925469510489, 1.4143413024221314], [0.6257047257907957, 0.3313089713994106, 1.3845218274856699], [0.613902676672381, 0.3376245167421365, 1.3555376100555918], [0.6022922233375773, 0.3435558171362999, 1.3273788786316545], [0.5908708630732223, 0.34911158622656124, 1.3000189023475366], [0.5796364066012643, 0.3543006260679922, 1.2734323579335558], [0.5685891427666308, 0.3591354585158367, 1.247603262337412], [0.5577244248972439, 0.3636198266556254, 1.222497780980065], [0.5470376817840281, 0.3677566985351186, 1.1980821655243636], [0.5365477011227113, 0.37158827960145374, 1.174408574998189], [0.5262518415954291, 0.375119254789794, 1.151450360605967], [0.5161117308644261, 0.3782898953496043, 1.1290471071478714], [0.506151235208352, 0.3811446245703845, 1.1072654466650917], [0.49643913783238963, 0.38381191064974546, 1.0863438981364582], [0.4869099973490608, 0.38616410247749555, 1.0660051430578972], [0.4775297524968833, 0.3881093941075793, 1.0460603253137397], [0.4684749785808654, 0.38992138052458586, 1.0271141081071267], [0.45974526665258986, 0.3915268548388567, 1.0091369350163977], [0.45133711875256377, 0.39268839845583026, 0.9918106442511835], [0.44360462088227315, 0.3936271347864685, 0.9758868450205826], [0.43694688490010786, 0.3945161021216901, 0.9624515377326197], [0.4315906996319098, 0.3950046644434845, 0.9519181120978144], [0.4280577796380111, 0.3950927169914573, 0.9449102789173061], [0.42662713526892415, 0.39523735049784253, 0.9418549215335764], [0.4263284481885641, 0.39539983660565586, 0.9411670544600845], [0.4261010039232765, 0.3954292866956077, 0.9407657792664712], [0.42572300397647544, 0.3954362409325253, 0.9401247690289407], [0.4254822282362977, 0.39531183877733356, 0.9397384841766858], [0.4254466207519215, 0.39530922766882426, 0.939649045403381], [0.42546940588572446, 0.3953262638504884, 0.9396702557447818], [0.4255639512761963, 0.3953617982565985, 0.9398113196777357], [0.4257198583352493, 0.3954116346724917, 0.9400546614968217], [0.42581302317993885, 0.3954067551524103, 0.9402653069304043], [0.4258526862271902, 0.39535395675278007, 0.9404027855577269], [0.4258714240641611, 0.395284462949466, 0.9404953011371568], [0.425892998601692, 0.3952400707410858, 0.9406011864423427], [0.4259409934742932, 0.3952421433588312, 0.9407357626952237], [0.42602463468104246, 0.39528491235281205, 0.940862658226654], [0.4260627358789707, 0.39528396093072804, 0.9409225413665789], [0.42605917275031174, 0.39524651805703004, 0.940933940865112], [0.4260409905508368, 0.3952063835017763, 0.9409455426049955], [0.4260143218690894, 0.3951704328543256, 0.9409706149125018], [0.42598504047492536, 0.3951382777685758, 0.9409968677473666], [0.42595361685841504, 0.39510373642457264, 0.9410066603631775], [0.4259176683861502, 0.39506249158531886, 0.9409949631469665], [0.4258735264798387, 0.39501312826939383, 0.9409692975433177], [0.4258160366423319, 0.39495457353098556, 0.9409429163760165], [0.4257263757082994, 0.39487008479951224, 0.9409124779135525], [0.4255673803997763, 0.3947237359440504, 0.9408577130251474], [0.42528467198193687, 0.3944625034470147, 0.9407470559057316], [0.4247993631905921, 0.39401166713552305, 0.9405422507817573], [0.4239920781342528, 0.3932607721085916, 0.940193657571234], [0.42267149025773365, 0.39203346864510247, 0.9396234505655598], [0.42052899514731706, 0.39004492442390837, 0.9387037645781037], [0.4170910329921099, 0.38685654318792917, 0.9372317661740734], [0.4116937582465378, 0.3818526534013186, 0.9349206407141326], [0.4034986308818913, 0.3742538331307116, 0.9314039539431204], [0.3915748458224322, 0.3631948519634849, 0.9262764741886143], [0.3750716594135053, 0.34788737415834053, 0.9191768559333721], [0.3534582917256142, 0.32783926350332976, 0.9098768363407757], [0.3274107131496261, 0.3036783682624417, 0.8986690284294172], [0.30408453946573677, 0.2820425956800543, 0.8886350519978191], [0.2880797380240309, 0.26719804652338947, 0.8817523154412776], [0.2777536023126187, 0.2576212632823299, 0.8773150265869698], [0.2714941800750227, 0.25181736637494, 0.8746309338916143], [0.26796455746462006, 0.2485463958292206, 0.8731252306655255], [0.2661489075306887, 0.2468649134704301, 0.8723559014706612], [0.2653321047138787, 0.24610848534969482, 0.8720105743455494], [0.26504501200525343, 0.24584261229239654, 0.871889837227163], [0.26499680536031, 0.24579806320220626, 0.8718705594723971], [0.26503028667104, 0.2458291480410418, 0.8718863536216314], [0.2651461179622153, 0.24593677737112146, 0.8719377869308926], [0.26532705049444316, 0.24610502136903495, 0.8720182829644928], [0.26546404993930506, 0.24623292450646928, 0.8720816510574851], [0.26552823985821244, 0.24629394823860584, 0.8721157447995378], [0.26555668912454156, 0.24632207841137932, 0.8721372782566016], [0.2655650612382585, 0.24633212864435414, 0.8721475540669055], [0.26556831597868324, 0.24633520464693826, 0.8721565732314354], [0.2655664071619714, 0.24633602392438264, 0.8721556421933657], [0.26556791706467286, 0.2463353499476911, 0.8721673128163995], [0.2655634097000559, 0.246336945844468, 0.8721585322545238], [0.2655695514269718, 0.24633486366368795, 0.872181468557377], [0.26555846164075636, 0.24633708149196326, 0.8721433990469964], [0.2655743647087028, 0.24632596868104287, 0.8721857625525193], [0.2655404231657049, 0.24632285954308838, 0.8720552538796724], [0.26557741583663935, 0.24628025952245067, 0.8721392193908999], [0.26548130566354944, 0.2462562709785318, 0.8717662012972478], [0.2655350148913803, 0.24611680810730002, 0.8718176131755311], [0.26521597756224985, 0.24569618264428061, 0.8702821434806106], [0.26463337078577814, 0.24413630980293802, 0.8667536454186975], [0.2623016973949709, 0.23933052770835972, 0.8540286121550217], [0.25558961503709754, 0.22456711313627545, 0.8168288786238507], [0.23621499150763098, 0.18450501347727138, 0.714197885466757], [0.19365227250474207, 0.10189335393450213, 0.5087787480888863], [0.14416508151782162, 0.02393096180894112, 0.3118825601635233], [0.12860573238610024, 0.004045930945053048, 0.2606316475394701], [0.12560082507923204, 0.0006467647375340353, 0.2517086873858034], [0.12509627016853242, 0.00010235023573430364, 0.2502707253156842], [0.12501526465877277, 1.6173958411031685e-5, 0.25004279053570444], [0.12500241433627954, 2.555884284393032e-6, 0.250006762190588], [0.12500038161702565, 4.0389755369405554e-7, 0.2500010686113911], [0.12500006030751598, 6.382476840976996e-8, 0.2500001688644371], [0.12500000952978313, 1.0085426613049774e-8, 0.2500000266835297], [0.12500000150583523, 1.5936285551910505e-9, 0.25000000421634455], [0.12500000023793517, 2.518071839296561e-10, 0.25000000066621897], [0.12500000003759457, 3.9786499787400154e-11, 0.2500000001052649], [0.1250000000059398, 6.286208811176184e-12, 0.2500000000166316], [0.1250000000009383, 9.930924953561664e-13, 0.25000000000262723], [0.1250000000001481, 1.5684124748503955e-13, 0.2500000000004148], [0.12500000000002326, 2.4719858524200256e-14, 0.2500000000000653], [0.12500000000000358, 3.853966538786583e-15, 0.25000000000001], [0.12500000000000053, 5.58211485729034e-16, 0.2500000000000014], [0.125, 6.665593845918698e-17, 0.2500000000000001], [0.125, 1.0806192927125272e-17, 0.25], [0.125, 0.0, 0.25], [0.125, 0.0, 0.25], [0.125, 0.0, 0.25], [0.125, 0.0, 0.25], [0.125, 0.0, 0.25], [0.125, -2.145979746602274e-18, 0.25], [0.125, -1.3744084525997318e-17, 0.25], [0.125, -1.478662382609276e-17, 0.25], [0.125, -1.388581318789706e-17, 0.25]], 0.2)Exact Solution
The exact Riemann solution for the Sod problem consists of five constant states separated by a rarefaction, contact, and shock. The pre-computed star-region values are:
function sod_exact(x, t; x0 = 0.5, gamma = 1.4)
rhoL, vL, PL = 1.0, 0.0, 1.0
rhoR, vR, PR = 0.125, 0.0, 0.1
cL = sqrt(gamma * PL / rhoL)
# Pre-computed star-region values
P_star = 0.30313017805064707
v_star = 0.92745262004895057
rho_star_L = 0.42631942817849544
rho_star_R = 0.26557371170530708
c_star_L = sqrt(gamma * P_star / rho_star_L)
# Shock speed
S_shock = v_star + 1.0 / (gamma * rho_star_R) *
(P_star - PR) / (v_star - vR + 1.0e-30)
xi = (x - x0) / t
if xi < -cL
return rhoL, vL, PL
elseif xi < v_star - c_star_L
# Inside rarefaction fan
v_fan = 2.0 / (gamma + 1) * (cL + xi)
c_fan = cL - 0.5 * (gamma - 1) * v_fan
rho_fan = rhoL * (c_fan / cL)^(2.0 / (gamma - 1))
P_fan = PL * (rho_fan / rhoL)^gamma
return rho_fan, v_fan, P_fan
elseif xi < v_star
return rho_star_L, v_star, P_star
elseif xi < S_shock
return rho_star_R, v_star, P_star
else
return rhoR, vR, PR
end
endsod_exact (generic function with 1 method)Visualisation
We now compare the numerical solutions against the exact solution.
using CairoMakie
x_exact = range(0.0, 1.0, length = 1000)
rho_exact = [sod_exact(xi, 0.2)[1] for xi in x_exact]
v_exact = [sod_exact(xi, 0.2)[2] for xi in x_exact]
P_exact = [sod_exact(xi, 0.2)[3] for xi in x_exact]
# Extract primitive variables from the numerical solutions
rho_hllc = [conserved_to_primitive(law, U_hllc[i])[1] for i in eachindex(U_hllc)]
v_hllc = [conserved_to_primitive(law, U_hllc[i])[2] for i in eachindex(U_hllc)]
P_hllc = [conserved_to_primitive(law, U_hllc[i])[3] for i in eachindex(U_hllc)]
rho_hll = [conserved_to_primitive(law, U_hll[i])[1] for i in eachindex(U_hll)]
v_hll = [conserved_to_primitive(law, U_hll[i])[2] for i in eachindex(U_hll)]
P_hll = [conserved_to_primitive(law, U_hll[i])[3] for i in eachindex(U_hll)]
fig = Figure(fontsize = 24, size = (1200, 400))
ax1 = Axis(fig[1, 1], xlabel = "x", ylabel = L"\rho", title = "Density")
ax2 = Axis(fig[1, 2], xlabel = "x", ylabel = "v", title = "Velocity")
ax3 = Axis(fig[1, 3], xlabel = "x", ylabel = "P", title = "Pressure")
lines!(ax1, x_exact, rho_exact, color = :black, linewidth = 2, label = "Exact")
scatter!(ax1, x_hllc, rho_hllc, color = :blue, markersize = 4, label = "HLLC")
scatter!(ax1, x_hll, rho_hll, color = :red, markersize = 4, label = "HLL")
lines!(ax2, x_exact, v_exact, color = :black, linewidth = 2)
scatter!(ax2, x_hllc, v_hllc, color = :blue, markersize = 4)
scatter!(ax2, x_hll, v_hll, color = :red, markersize = 4)
lines!(ax3, x_exact, P_exact, color = :black, linewidth = 2)
scatter!(ax3, x_hllc, P_hllc, color = :blue, markersize = 4)
scatter!(ax3, x_hll, P_hll, color = :red, markersize = 4)
axislegend(ax1, position = :cb)
resize_to_layout!(fig)
fig
Notice that the HLLC solver (blue) resolves the contact discontinuity much more sharply than HLL (red), which smears it out. Both solvers capture the shock and rarefaction well at this resolution ($N = 200$).
Convergence
We can verify that the error decreases with resolution:
function compute_l1_error(N)
m = StructuredMesh1D(0.0, 1.0, N)
p = HyperbolicProblem(
law, m, HLLCSolver(), CellCenteredMUSCL(MinmodLimiter()),
bc_left, bc_right, ic; final_time = 0.2, cfl = 0.5
)
xx, UU, _ = solve_hyperbolic(p)
err = 0.0
for i in eachindex(xx)
rho_num = conserved_to_primitive(law, UU[i])[1]
rho_ex = sod_exact(xx[i], 0.2)[1]
err += abs(rho_num - rho_ex)
end
return err / N
end
err_100 = compute_l1_error(100)
err_200 = compute_l1_error(200)
err_400 = compute_l1_error(400)trueJust the code
An uncommented version of this example is given below. You can view the source code for this file here.
using FiniteVolumeMethod
using StaticArrays
gamma = 1.4
eos = IdealGasEOS(gamma)
law = EulerEquations{1}(eos)
wL = SVector(1.0, 0.0, 1.0)
wR = SVector(0.125, 0.0, 0.1)
N = 200
mesh = StructuredMesh1D(0.0, 1.0, N)
bc_left = DirichletHyperbolicBC(wL)
bc_right = DirichletHyperbolicBC(wR)
ic(x) = x < 0.5 ? wL : wR
prob_hllc = HyperbolicProblem(
law, mesh, HLLCSolver(), CellCenteredMUSCL(MinmodLimiter()),
bc_left, bc_right, ic;
final_time = 0.2, cfl = 0.5
)
x_hllc, U_hllc, t_hllc = solve_hyperbolic(prob_hllc)
prob_hll = HyperbolicProblem(
law, mesh, HLLSolver(), CellCenteredMUSCL(MinmodLimiter()),
bc_left, bc_right, ic;
final_time = 0.2, cfl = 0.5
)
x_hll, U_hll, t_hll = solve_hyperbolic(prob_hll)
function sod_exact(x, t; x0 = 0.5, gamma = 1.4)
rhoL, vL, PL = 1.0, 0.0, 1.0
rhoR, vR, PR = 0.125, 0.0, 0.1
cL = sqrt(gamma * PL / rhoL)
# Pre-computed star-region values
P_star = 0.30313017805064707
v_star = 0.92745262004895057
rho_star_L = 0.42631942817849544
rho_star_R = 0.26557371170530708
c_star_L = sqrt(gamma * P_star / rho_star_L)
# Shock speed
S_shock = v_star + 1.0 / (gamma * rho_star_R) *
(P_star - PR) / (v_star - vR + 1.0e-30)
xi = (x - x0) / t
if xi < -cL
return rhoL, vL, PL
elseif xi < v_star - c_star_L
# Inside rarefaction fan
v_fan = 2.0 / (gamma + 1) * (cL + xi)
c_fan = cL - 0.5 * (gamma - 1) * v_fan
rho_fan = rhoL * (c_fan / cL)^(2.0 / (gamma - 1))
P_fan = PL * (rho_fan / rhoL)^gamma
return rho_fan, v_fan, P_fan
elseif xi < v_star
return rho_star_L, v_star, P_star
elseif xi < S_shock
return rho_star_R, v_star, P_star
else
return rhoR, vR, PR
end
end
using CairoMakie
x_exact = range(0.0, 1.0, length = 1000)
rho_exact = [sod_exact(xi, 0.2)[1] for xi in x_exact]
v_exact = [sod_exact(xi, 0.2)[2] for xi in x_exact]
P_exact = [sod_exact(xi, 0.2)[3] for xi in x_exact]
# Extract primitive variables from the numerical solutions
rho_hllc = [conserved_to_primitive(law, U_hllc[i])[1] for i in eachindex(U_hllc)]
v_hllc = [conserved_to_primitive(law, U_hllc[i])[2] for i in eachindex(U_hllc)]
P_hllc = [conserved_to_primitive(law, U_hllc[i])[3] for i in eachindex(U_hllc)]
rho_hll = [conserved_to_primitive(law, U_hll[i])[1] for i in eachindex(U_hll)]
v_hll = [conserved_to_primitive(law, U_hll[i])[2] for i in eachindex(U_hll)]
P_hll = [conserved_to_primitive(law, U_hll[i])[3] for i in eachindex(U_hll)]
fig = Figure(fontsize = 24, size = (1200, 400))
ax1 = Axis(fig[1, 1], xlabel = "x", ylabel = L"\rho", title = "Density")
ax2 = Axis(fig[1, 2], xlabel = "x", ylabel = "v", title = "Velocity")
ax3 = Axis(fig[1, 3], xlabel = "x", ylabel = "P", title = "Pressure")
lines!(ax1, x_exact, rho_exact, color = :black, linewidth = 2, label = "Exact")
scatter!(ax1, x_hllc, rho_hllc, color = :blue, markersize = 4, label = "HLLC")
scatter!(ax1, x_hll, rho_hll, color = :red, markersize = 4, label = "HLL")
lines!(ax2, x_exact, v_exact, color = :black, linewidth = 2)
scatter!(ax2, x_hllc, v_hllc, color = :blue, markersize = 4)
scatter!(ax2, x_hll, v_hll, color = :red, markersize = 4)
lines!(ax3, x_exact, P_exact, color = :black, linewidth = 2)
scatter!(ax3, x_hllc, P_hllc, color = :blue, markersize = 4)
scatter!(ax3, x_hll, P_hll, color = :red, markersize = 4)
axislegend(ax1, position = :cb)
resize_to_layout!(fig)
fig
function compute_l1_error(N)
m = StructuredMesh1D(0.0, 1.0, N)
p = HyperbolicProblem(
law, m, HLLCSolver(), CellCenteredMUSCL(MinmodLimiter()),
bc_left, bc_right, ic; final_time = 0.2, cfl = 0.5
)
xx, UU, _ = solve_hyperbolic(p)
err = 0.0
for i in eachindex(xx)
rho_num = conserved_to_primitive(law, UU[i])[1]
rho_ex = sod_exact(xx[i], 0.2)[1]
err += abs(rho_num - rho_ex)
end
return err / N
end
err_100 = compute_l1_error(100)
err_200 = compute_l1_error(200)
err_400 = compute_l1_error(400)This page was generated using Literate.jl.