Skip to content
Snippets Groups Projects
Commit f9f6f793 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

dg : test 1D

parent 416425c0
No related branches found
No related tags found
No related merge requests found
...@@ -3,6 +3,7 @@ ...@@ -3,6 +3,7 @@
--]] --]]
x = os.clock()
function initial_condition( xyz , f ) function initial_condition( xyz , f )
for i=0,xyz:size1()-1 do for i=0,xyz:size1()-1 do
f:set (i, 0, 3.14) f:set (i, 0, 3.14)
...@@ -12,22 +13,35 @@ end ...@@ -12,22 +13,35 @@ end
PI = 3.14159 PI = 3.14159
T= 0.05 T= 0.05
t=0
Amax = 3.15 Amax = 3.15
function inlet( FCT ) function inlet( FCT )
FCT:set(0,0,Amax*math.sin(2*PI*t/T)*math.step(T/2-t)) FCT:set(0,0,Amax*math.sin(2*PI*t/T))
FCT:set(0,1, 0) FCT:set(0,1, 0)
end end
--[[
function pressureLawShallowWater (solution, bathymetry, f)
for i=0,f:size1()-1 do
f:set (i, 0, bathymetry:get(i,0))
end
end
--]]
xyz = functionCoordinates.get() xyz = functionCoordinates.get()
--[[ --[[
Example of a lua program driving the DG code Example of a lua program driving the DG code
--]] --]]
order=2
options = gmshOptions()
model = GModel() model = GModel()
model:load ('aorta.msh') model:load ('aorta.geo')
order=3 options:numberSet('Mesh',-1,'ElementOrder',order)
model:mesh(1)
model:save('aorta.msh')
dimension=1 dimension=1
-- boundary condition -- boundary condition
...@@ -35,7 +49,9 @@ inlet_bc = functionLua(2, 'inlet') ...@@ -35,7 +49,9 @@ inlet_bc = functionLua(2, 'inlet')
law = dgConservationLawShallowWater1d() law = dgConservationLawShallowWater1d()
law:addBoundaryCondition('Inlet',law:newOutsideValueBoundary(inlet_bc)) law:addBoundaryCondition('Inlet',law:newOutsideValueBoundary(inlet_bc))
law:addBoundaryCondition('Outlet',law:newBoundaryWall()) law:addBoundaryCondition('Outlet',law:newBoundaryWall())
law:setBathymetry(functionConstant({0})) bathymetry = functionConstant({0})
law:setBathymetry(bathymetry)
--law:setPressureLaw(functionLua(1,'pressureLaw',{functionSolution.get(), bathymetry}))
law:setLinearDissipation(functionConstant({0})) law:setLinearDissipation(functionConstant({0}))
groups = dgGroupCollection(model, dimension, order) groups = dgGroupCollection(model, dimension, order)
...@@ -44,6 +60,7 @@ groups:buildGroupsOfInterfaces() ...@@ -44,6 +60,7 @@ groups:buildGroupsOfInterfaces()
rk=dgRungeKutta() rk=dgRungeKutta()
limiter = dgSlopeLimiter(law) limiter = dgSlopeLimiter(law)
print'*** print projetion ***'
-- build solution vector -- build solution vector
FS = functionLua(2, 'initial_condition', {xyz}) FS = functionLua(2, 'initial_condition', {xyz})
solution = dgDofContainer(groups, law:getNbFields()) solution = dgDofContainer(groups, law:getNbFields())
...@@ -57,6 +74,7 @@ print'*** solve ***' ...@@ -57,6 +74,7 @@ print'*** solve ***'
CFL = 0.2; CFL = 0.2;
for i=1,9000 do for i=1,9000 do
dt = CFL * rk:computeInvSpectralRadius(law,solution); dt = CFL * rk:computeInvSpectralRadius(law,solution);
t = t+dt
norm = rk:iterate44(law,dt,solution) norm = rk:iterate44(law,dt,solution)
if (i % 100 == 0) then if (i % 100 == 0) then
print('|ITER|',i,'|NORM|',norm,'|DT|',dt,'|CPU|',os.clock() - x) print('|ITER|',i,'|NORM|',norm,'|DT|',dt,'|CPU|',os.clock() - x)
......
model = GModel() model = GModel()
model:load ('stommel_square.msh') model:load ('stommel_square.msh')
order=1 order = 2
dimension = 2 dimension = 2
CFunctions =[[ CFunctions =[[
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment