From 304facceb63887b6d7d5af3e43a474912abf07a0 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Fri, 19 Mar 2010 11:45:12 +0000
Subject: [PATCH]
---
Solver/TESTCASES/Aorta.lua | 92 ++++++++++++++++++++++++-----------
Solver/TESTCASES/DamBreak.lua | 49 ++++++++++++++-----
Solver/TESTCASES/aorta.geo | 6 +--
3 files changed, 105 insertions(+), 42 deletions(-)
diff --git a/Solver/TESTCASES/Aorta.lua b/Solver/TESTCASES/Aorta.lua
index 44aacd8bd9..d1fc49859a 100644
--- a/Solver/TESTCASES/Aorta.lua
+++ b/Solver/TESTCASES/Aorta.lua
@@ -1,68 +1,104 @@
--[[
Function for initial conditions
+ CGS UNITS
--]]
+t=0
+Ainit = 3.0
+T= 0.05
+Amax = 0.01
+p0 = 0.0
+rho = 1.06
+--c = 500
+--beta = 2*rho*math.pow(c,2)/math.sqrt(A0)
+beta = 282166
+--KR = 8*math.pi*0.04/rho
+KR = 0.0
+
+--[[
+ Functions
+--]]
-x = os.clock()
function initial_condition( xyz , f )
for i=0,xyz:size1()-1 do
- f:set (i, 0, 3.14)
+ f:set (i, 0, 0)
f:set (i, 1, 0)
end
end
-PI = 3.14159
-T= 0.05
-t=0
-Amax = 3.15
-
function inlet( FCT )
- FCT:set(0,0,Amax*math.sin(2*PI*t/T)*atan(1000*(T/2-t))
+ FCT:set(0,0,Amax*math.sin(2*math.pi*t/T)*math.atan(1000*(T/2-t)))
+ --FCT:set(0,0, Amax)
FCT:set(0,1, 0)
end
---[[
-function pressureLawShallowWater (solution, bathymetry, f)
+function pressureFlux (solution, bathymetry, f)
for i=0,f:size1()-1 do
- f:set (i, 0, bathymetry:get(i,0))
+ A0 = bathymetry:get(i,0)
+ A = A0 + solution:get(i,0)
+ p = p0 + beta*(math.sqrt(A)-math.sqrt(A0))
+ f:set (i, 0, p)
end
end
---]]
-xyz = functionCoordinates.get()
+function celerityFunction (solution, bathymetry, f)
+ for i=0,f:size1()-1 do
+ A0 = bathymetry:get(i,0)
+ c = math.sqrt(beta/(2*rho))*math.pow(A0, 0.25)
+ f:set (i, 0, c)
+ end
+end
+
+
+function characteristicFunction (solution, celerity, f)
+ for i=0,f:size1()-1 do
+ W1 = sol:get(i,1) + 4*celerity:get(i,0)
+ W2 = sol:get(i,1) - 4*celerity:get(i,0)
+ f:set (i, 0, W1)
+ f:set (i, 1, W2)
+ end
+end
--[[
Example of a lua program driving the DG code
--]]
-order=2
-options = gmshOptions()
-model = GModel()
-model:load ('aorta.geo')
-options:numberSet('Mesh',-1,'ElementOrder',order)
-model:mesh(1)
-model:save('aorta.msh')
+order=1
dimension=1
+model = GModel()
+model:load('aorta.msh')
+--options = gmshOptions()
+--model:load ('aorta.geo')
+--options:numberSet('Mesh',-1,'ElementOrder',order)
+--model:mesh(dimension)
+--model:save('aorta.msh')
+
+bathymetry = functionConstant({Ainit})
+dissipation = functionConstant({KR})
+pressure = functionLua(1,'pressureFlux', {functionSolution.get(), bathymetry})
+pwv = functionLua(1,'celerityFunction', {functionSolution.get(), bathymetry})
+--chars = functionLua(1,'characteristicFunction', {functionSolution.get(), pwv})
+
-- boundary condition
inlet_bc = functionLua(2, 'inlet')
law = dgConservationLawShallowWater1d()
law:addBoundaryCondition('Inlet',law:newOutsideValueBoundary(inlet_bc))
law:addBoundaryCondition('Outlet',law:newBoundaryWall())
-bathymetry = functionConstant({0})
+law:setLinearDissipation(dissipation)
law:setBathymetry(bathymetry)
---law:setPressureLaw(functionLua(1,'pressureLaw',{functionSolution.get(), bathymetry}))
-law:setLinearDissipation(functionConstant({0}))
+law:setPressureFlux(pressure)
+law:setCelerity(pwv)
+--law:setCharacteristicFunction(chars)
groups = dgGroupCollection(model, dimension, order)
groups:buildGroupsOfInterfaces()
rk=dgRungeKutta()
-limiter = dgSlopeLimiter(law)
print'*** print projetion ***'
-- build solution vector
-FS = functionLua(2, 'initial_condition', {xyz})
+FS = functionLua(2, 'initial_condition', {functionCoordinates.get()})
solution = dgDofContainer(groups, law:getNbFields())
solution:L2Projection(FS)
@@ -70,14 +106,14 @@ print'*** print initial sol ***'
solution:exportMsh('output/init')
print'*** solve ***'
---dt = 0.00001;
CFL = 0.2;
-for i=1,9000 do
+x = os.clock()
+for i=1,10000 do
dt = CFL * rk:computeInvSpectralRadius(law,solution);
t = t+dt
norm = rk:iterate44(law,dt,solution)
if (i % 100 == 0) then
- print('|ITER|',i,'|NORM|',norm,'|DT|',dt,'|CPU|',os.clock() - x)
+ print('|ITER|',i,'|NORM|',norm,'|DT|',dt, '|t|', t, '|CPU|',os.clock() - x)
end
if (i % 100 == 0) then
solution:exportMsh(string.format('output/solution-%06d', i))
diff --git a/Solver/TESTCASES/DamBreak.lua b/Solver/TESTCASES/DamBreak.lua
index 52bc169308..cf68af5161 100644
--- a/Solver/TESTCASES/DamBreak.lua
+++ b/Solver/TESTCASES/DamBreak.lua
@@ -7,8 +7,6 @@ g = 9.81;
function initial_condition( xyz , f )
for i=0,xyz:size1()-1 do
x = xyz:get(i,0)
- y = xyz:get(i,1)
- z = xyz:get(i,2)
if (x<0.0) then
f:set (i, 0, 40)
f:set (i, 1, 0)
@@ -19,19 +17,47 @@ function initial_condition( xyz , f )
end
end
-function pressureFlux (solution, f)
+function pressureFlux (solution, f)
for i=0,f:size1()-1 do
p = g*solution:get(i,0)
f:set (i, 0, p)
end
end
-function celerity (solution, bathymetry, f)
+function celerityFunction (solution, bathymetry, f)
for i=0,f:size1()-1 do
c = math.sqrt(g*(bathymetry:get(i,0)+solution:get(i,0)))
- f:set (i, 0, c)
+ f:set(i, 0, c)
end
end
+
+--CFunctions =[[
+--void pressureFlux ( fullMatrix<double> &f, fullMatrix<double> &sol) {
+-- for (size_t i = 0; i< f.size1(); i++) {
+-- double p=9.81*sol(i,0);
+-- f.set(i,0,p);
+-- f.set(i,1,0);
+-- }
+--}
+--void celerityFunction (fullMatrix<double> &f , fullMatrix<double> &sol, fullMatrix<double> &bathymetry) {
+-- for (size_t i = 0; i< sol.size1(); i++) {
+-- double c =sqrt(9.81*(bathymetry(i,0)+sol(i,0)));
+-- sol.set(i,0,c);
+-- }
+--}]]
+--if (Msg.getCommRank() == 0 ) then
+-- cfile = io.popen("g++-4.2 -O3 -pipe -m32 -shared -o tmp.dylib -I ../../Numeric -I../../Common -I../../bin/Common -x c++ - ","w");
+-- cfile:write("#include\"fullMatrix.h\"\nextern \"C\" {")
+-- cfile:write(CFunctions)
+-- cfile:write("}")
+-- cfile:close()
+--end
+--Msg.barrier()
+
+
+bathymetry = functionConstant({0})
+dissipation = functionConstant({0})
+
--[[
Example of a lua program driving the DG code
--]]
@@ -41,8 +67,6 @@ model:load ('edge.msh')
order=1
dimension=1
-bathymetry = functionConstant({0})
-dissipation = functionConstant({0})
-- boundary condition
law = dgConservationLawShallowWater1d()
@@ -51,7 +75,9 @@ law:addBoundaryCondition('Right',law:newBoundaryWall())
law:setLinearDissipation(dissipation)
law:setBathymetry(bathymetry)
law:setPressureFlux(functionLua(1,'pressureFlux',{functionSolution.get()}))
-law:setCelerity(functionLua(1,'celerity',{functionSolution.get(), bathymetry}))
+law:setCelerity(functionLua(1,'celerityFunction',{functionSolution.get(), bathymetry}))
+--law:setPressureFlux(functionC("tmp.dylib","pressureFlux",1,{functionSolution.get()}))
+--law:setPressureFlux(functionC("tmp.dylib","celerityFunction",1,{functionSolution.get(), bathymetry}))
groups = dgGroupCollection(model, dimension, order)
groups:buildGroupsOfInterfaces()
@@ -61,8 +87,7 @@ limiter = dgSlopeLimiter(law)
rk:setLimiter(limiter)
-- build solution vector
-xyz = functionCoordinates.get()
-FS = functionLua(2, 'initial_condition', {xyz})
+FS = functionLua(2, 'initial_condition', {functionCoordinates.get()})
solution = dgDofContainer(groups, law:getNbFields())
solution:L2Projection(FS)
@@ -72,7 +97,9 @@ limiter:apply(solution)
solution:exportMsh('output/init_limit')
print'*** solve ***'
---dt = 0.00001;
+
+
+
CFL = 0.2;
for i=1,15000 do
dt = CFL * rk:computeInvSpectralRadius(law,solution);
diff --git a/Solver/TESTCASES/aorta.geo b/Solver/TESTCASES/aorta.geo
index 70ea0e186a..bb78005fd6 100644
--- a/Solver/TESTCASES/aorta.geo
+++ b/Solver/TESTCASES/aorta.geo
@@ -1,6 +1,6 @@
-
-Point(1) = {0, 0, 0, 0.5};
-Point(2) = {100, 0, 0, 0.5};
+lc = 0.5*50;
+Point(1) = {0, 0, 0, lc};
+Point(2) = {100, 0, 0, lc};
Line(1) = {1, 2};
Physical Point("Inlet") = {1};
--
GitLab