No commit message

No commit message

Solver/TESTCASES/Aorta.lua

 --[[ 
      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
+
+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
---]]
 
-xyz = functionCoordinates.get()
+
+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)) 

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)
@@ -25,13 +23,41 @@ function pressureFlux (solution, f)
      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)
   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);  

Solver/TESTCASES/aorta.geo

-
-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};