Skip to content
Snippets Groups Projects
Commit bb5dea91 authored by Emilie Marchandise's avatar Emilie Marchandise
Browse files

No commit message

No commit message
parent e1b2b8c1
No related branches found
No related tags found
No related merge requests found
...@@ -4,7 +4,7 @@ model:load ('edge.geo') ...@@ -4,7 +4,7 @@ model:load ('edge.geo')
-- model:save ('edge.msh') -- model:save ('edge.msh')
model:load ('edge.msh') model:load ('edge.msh')
dg = dgSystemOfEquations (model) dg = dgSystemOfEquations (model)
dg:setOrder(0) dg:setOrder(1)
-- conservation law -- conservation law
-- advection speed -- advection speed
...@@ -54,7 +54,7 @@ dg:exportSolution('output/Adv1D-00000') ...@@ -54,7 +54,7 @@ dg:exportSolution('output/Adv1D-00000')
-- main loop -- main loop
n = 5 n = 5
for i=1,100*n do for i=1,100*n do
norm = dg:RK44(0.03) norm = dg:RK44_limiter(0.03)
if (i % n == 0) then if (i % n == 0) then
print('iter',i,norm) print('iter',i,norm)
dg:exportSolution(string.format("output/Adv1D-%05d", i)) dg:exportSolution(string.format("output/Adv1D-%05d", i))
......
MACH = 3; MACH = 1.0;
RHO = 1.0; GAMMA = 1.4;
PRES = 1./(MACH*RHO*RHO*1.4*1.4) U = 3.0
V = 1.0 V = 0.0
SOUND = V/MACH RHO = 1.4;
PRES = RHO*U*U/(GAMMA*MACH*MACH)
--PRES = 1;
--PRES = ./(MACH*RHO*RHO*GAMMA*GAMMA)
SOUND = math.sqrt(U*U+V*V)/MACH
--[[ --[[
Function for initial conditions Function for initial conditions
...@@ -10,9 +16,9 @@ SOUND = V/MACH ...@@ -10,9 +16,9 @@ SOUND = V/MACH
function free_stream( XYZ, FCT ) function free_stream( XYZ, FCT )
for i=0,XYZ:size1()-1 do for i=0,XYZ:size1()-1 do
FCT:set(i,0,RHO) FCT:set(i,0,RHO)
FCT:set(i,1,RHO*V) FCT:set(i,1,RHO*U)
FCT:set(i,2,0.0) FCT:set(i,2,RHO*V)
FCT:set(i,3, 0.5*RHO*V*V+PRES/0.4) FCT:set(i,3, 0.5*RHO*(U*U+V*V)+PRES/(RHO*GAMMA-1))
end end
end end
...@@ -25,18 +31,19 @@ print'*** Loading the mesh and the model ***' ...@@ -25,18 +31,19 @@ print'*** Loading the mesh and the model ***'
myModel = GModel () myModel = GModel ()
myModel:load ('step.geo') myModel:load ('step.geo')
myModel:load ('step.msh') myModel:load ('step.msh')
print'*** Create a dg solver ***' print'*** Create a dg solver ***'
DG = dgSystemOfEquations (myModel) DG = dgSystemOfEquations (myModel)
DG:setOrder(order) DG:setOrder(order)
FS = functionLua(4, 'free_stream', {'XYZ'}):getName()
law=dgPerfectGasLaw2d() law=dgPerfectGasLaw2d()
DG:setConservationLaw(law) DG:setConservationLaw(law)
law:addBoundaryCondition('Walls',law:newWallBoundary())
FS = functionLua(4, 'free_stream', {'XYZ'}):getName()
law:addBoundaryCondition('Walls',law:newWallBoundary())
law:addBoundaryCondition('LeftRight',law:newOutsideValueBoundary(FS)) law:addBoundaryCondition('LeftRight',law:newOutsideValueBoundary(FS))
DG:setup()
DG:setup()
print'*** setting the initial solution ***' print'*** setting the initial solution ***'
...@@ -53,10 +60,10 @@ LC = 0.1*.1 ...@@ -53,10 +60,10 @@ LC = 0.1*.1
dt = .3*LC/(SOUND+V); dt = .3*LC/(SOUND+V);
print('DT=',dt) print('DT=',dt)
for i=1,10000 do for i=1,10 do
norm = DG:RK44_limiter(dt) norm = DG:RK44_limiter(0.1*dt)
print('*** ITER ***',i,norm) print('*** ITER ***',i,norm)
if (i % 100 == 0) then if (i % 1 == 0) then
DG:exportSolution(string.format("output/solution-%06d", i)) DG:exportSolution(string.format("output/solution-%06d", i))
end end
end end
......
Point(1) = {0, 0, 0, .1}; Point(1) = {0, 0, 0, .1};
Point(2) = {.5, 0, 0, .1}; Point(2) = {.6, 0, 0, .1};
Point(3) = {.5, .25, 0, .1}; Point(3) = {.6, .2, 0, .1};
Point(4) = {3, .25, 0, .1}; Point(4) = {2.4, .2, 0, .1};
Point(5) = {3, 1, 0, .1}; Point(5) = {2.4, 1, 0, .1};
Point(6) = {0, 1, 0, .1}; Point(6) = {0, 1, 0, .1};
Line(1) = {6, 5}; Line(1) = {6, 5};
Line(2) = {5, 4}; Line(2) = {5, 4};
...@@ -14,5 +14,4 @@ Line Loop(7) = {2, 3, 4, 5, 6, 1}; ...@@ -14,5 +14,4 @@ Line Loop(7) = {2, 3, 4, 5, 6, 1};
Plane Surface(8) = {7}; Plane Surface(8) = {7};
Physical Line("Walls") = {1, 3, 4, 5}; Physical Line("Walls") = {1, 3, 4, 5};
Physical Line("LeftRight") = {2, 6}; Physical Line("LeftRight") = {2, 6};
//Physical Line("LeftRight") = {1,2,3,4,5,6};
Physical Surface("Body") = {8}; Physical Surface("Body") = {8};
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment