Skip to content

Commit

Permalink
PDE extension
Browse files Browse the repository at this point in the history 
experimental code

git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk/doc@22287 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
Jan Silar committed Sep 12, 2014
1 parent 060ef56 commit 3c23445
Show file tree
Hide file tree
Showing 7 changed files with 315 additions and 86 deletions.
25 changes: 25 additions & 0 deletions PDEInModelica/julia/advTest.jl
View file
@@ -0,0 +1,25 @@
include("advection.jl")
tEnd = 0.5
(X,U) = simulate(tEnd)
nX = size(X,1)
function analyticSol(x,t)
if x >= c*t
initFun(1, x - c*t)
else
lBCFun(t - x/c)
end
end
analyticU = Array(Float64,nX)
for i in 1:nX
analyticU[i] = analyticSol(X[i],tEnd)
end


using Gadfly
#using Winston
#plot(x = X, y = U[1,:])


#D = vec(U[1,:]) - analyticU
#plot(x = X, y = D)
plot(x = X, y = analyticU)
118 changes: 32 additions & 86 deletions PDEInModelica/julia/advection.jl
View file
@@ -1,106 +1,52 @@
#model pder(u,t) + 1*pder(u,x) = 0
#model pder(u,t) + c*pder(u,x) = 0
include("incl.jl")

# - model static data
nU = 1 #number of state fields


#numerics setup
nX = 1000 #number of nodes

# - model values
L = 1.0 #length of domain
tEnd = 0.5

c = 1.0

#model functions:

function initModelFun(x,u)
for i = 1:size(u,2)
u[1,i] = 0.0
function initFun(i,x)
if i == 1
0.0
else
error("wrong state number in advection")
end
end

function evalBoundariesFun(u,t)
nX = size(u,2)
u[1] = sin(2.0*pi*t)
u[nX] = extrapolate3(x[nX],x[nX-1],x[nX-2],x[nX-3],u[nX-1],u[nX-2],u[nX-3])
end

function maxEigValFun()
1
end

function utFun(x,u,ux,t)
-1*ux
end


#numerics variables
x = linspace(0.0,L,nX) #coordinate
u = Array(Float64,nU,nX) #state fields u[nVar, nNode]
ux = Array(Float64,nU,nX) #state fields space derivatives ux[nVar, nNode]
ut = Array(Float64,nU,nX) #state fields time derivatives ut[nVar, nNode]
t = 0.0 #time
dx = L/(nX-1) #space step
#dt #time step
cfl = 0.5

#numerics functions
function extrapolate3(x,x1,x2,x3,u1,u2,u3)
#lagrange polinomial extrapolation
# q1 = u1 / ( (x1-x2)*(x1-x3) )
# q2 = u2 / ( (x2-x1)*(x2-x3) )
# q3 = u3 / ( (x3-x1)*(x3-x2) )
# x^2*(q1+q2+q3) + x*( q1*(x2+x3) + q2*(x1+x3) + q3*(x1+x2) ) + q1*x2*x3 + q2*x1*x3 + q3*x1*x2
u1*(x-x2)*(x-x3)/((x1-x2)*(x1-x3)) +
u2*((x-x1)*(x-x3))/((x2-x1)*(x2-x3)) +
u3*((x-x1)*(x-x2))/((x3-x1)*(x3-x2))
end

function updateUxLW(x,u,ux)
for i = 2:(size(x,1)-1)
ux[i] = (u[i+1] - u[i-1])/(x[i+1] - x[i-1])
function lBCFun(t)
if 0 < t < 0.1
sin(10.0*2.0*pi*t)
else
0
end
end

function updateUt(x,u,ux,ut,t)
for i = 2:size(x,1)-1
ut[i] = utFun(x[i],u[i],ux[i],t)
function BCFun(nState,side,t,X,U)
if nState == 1
if side == left lBCFun(t)
elseif side == right extrapolate(X,U,nState,side)
end
else
error("wrong state index in advection")
end
end

function updateULW(u,ut,dt)
function UUpdate(i)
(u[i+1] + u[i-1])/2 + dt*ut[i]
end
u1new = UUpdate(2)
unew = UUpdate(3)
for i = 4:size(u,2)-1
u[i-2] = u1new
u1new = unew
unew = UUpdate(i)
end
u[size(u,2)-2] = u1new
u[size(u,2)-1] = unew
function maxEigValFun()
c
end






#computation
using Gadfly
initModelFun(x,u)
while t < tEnd
dt = cfl*maxEigValFun()*dx
evalBoundariesFun(u,t)
updateUxLW(x,u,ux)
updateUt(x,u,ux,ut,t)
updateULW(u,ut,dt)
t = t + dt
# println(t)
end
plot(x = x, y = u[1,:])



function utFun(x,u,ux,t,iState)
if iState == 1
-c*ux
else
error("wrong variable index in utFun()")
end
end

include("solver.jl")
110 changes: 110 additions & 0 deletions PDEInModelica/julia/advectionSolver.jl
View file
@@ -0,0 +1,110 @@
#model pder(u,t) + 1*pder(u,x) = 0
# - model static data
nU = 1 #number of state fields


#numerics setup
nX = 100 #number of nodes

# - model values
L = 1.0 #length of domain
tEnd = 0.5

#model functions:

function initModelFun(x,u)
for i = 1:size(u,2)
u[1,i] = 0.0
end
end

function evalBoundariesFun(u,t)
nX = size(u,2)
u[1] = if 0 < t < 0.1
sin(10.0*2.0*pi*t)
else
0
end
u[nX] = extrapolate3(x[nX],x[nX-1],x[nX-2],x[nX-3],u[nX-1],u[nX-2],u[nX-3])
end

function maxEigValFun()
1
end

function utFun(x,u,ux,t)
-1*ux
end


#numerics variables
x = linspace(0.0,L,nX) #coordinate
u = Array(Float64,nU,nX) #state fields u[nVar, nNode]
ux = Array(Float64,nU,nX) #state fields space derivatives ux[nVar, nNode]
ut = Array(Float64,nU,nX) #state fields time derivatives ut[nVar, nNode]
t = 0.0 #time
dx = L/(nX-1) #space step
#dt #time step
cfl = 0.1

#numerics functions
function extrapolate3(x,x1,x2,x3,u1,u2,u3)
#lagrange polinomial extrapolation
# q1 = u1 / ( (x1-x2)*(x1-x3) )
# q2 = u2 / ( (x2-x1)*(x2-x3) )
# q3 = u3 / ( (x3-x1)*(x3-x2) )
# x^2*(q1+q2+q3) + x*( q1*(x2+x3) + q2*(x1+x3) + q3*(x1+x2) ) + q1*x2*x3 + q2*x1*x3 + q3*x1*x2
u1*(x-x2)*(x-x3)/((x1-x2)*(x1-x3)) +
u2*((x-x1)*(x-x3))/((x2-x1)*(x2-x3)) +
u3*((x-x1)*(x-x2))/((x3-x1)*(x3-x2))
end

function updateUxLW(x,u,ux)
for i = 2:(size(x,1)-1)
ux[i] = (u[i+1] - u[i-1])/(x[i+1] - x[i-1])
end
end

function updateUt(x,u,ux,ut,t)
for i = 2:size(x,1)-1
ut[i] = utFun(x[i],u[i],ux[i],t)
end
end

function updateULW(u,ut,dt)
function UUpdate(i)
(u[i+1] + u[i-1])/2 + dt*ut[i]
end
u1new = UUpdate(2)
unew = UUpdate(3)
for i = 4:size(u,2)-1
u[i-2] = u1new
u1new = unew
unew = UUpdate(i)
end
u[size(u,2)-2] = u1new
u[size(u,2)-1] = unew
end






#computation
using Gadfly
initModelFun(x,u)
while t < tEnd
dt = cfl*maxEigValFun()*dx
evalBoundariesFun(u,t)
updateUxLW(x,u,ux)
updateUt(x,u,ux,ut,t)
updateULW(u,ut,dt)
t = t + dt
# println(t)
end
plot(x = x, y = u[1,:])




21 changes: 21 additions & 0 deletions PDEInModelica/julia/incl.jl
View file
@@ -0,0 +1,21 @@
const left = false
const right = true

function extrapolate(X,U,nState,side)
if side == left
extrapolate3(X[1],X[2],X[3],X[4],U[nState,2],U[nState,3],U[nState,4])
elseif side == right
extrapolate3(X[nX],X[nX-1],X[nX-2],X[nX-3],U[nState,nX-1],U[nState,nX-2],U[nState,nX-3])
end
end

function extrapolate3(x,x1,x2,x3,u1,u2,u3)
#lagrange polinomial extrapolation
# q1 = u1 / ( (x1-x2)*(x1-x3) )
# q2 = u2 / ( (x2-x1)*(x2-x3) )
# q3 = u3 / ( (x3-x1)*(x3-x2) )
# x^2*(q1+q2+q3) + x*( q1*(x2+x3) + q2*(x1+x3) + q3*(x1+x2) ) + q1*x2*x3 + q2*x1*x3 + q3*x1*x2
u1*(x-x2)*(x-x3)/((x1-x2)*(x1-x3)) +
u2*((x-x1)*(x-x3))/((x2-x1)*(x2-x3)) +
u3*((x-x1)*(x-x2))/((x3-x1)*(x3-x2))
end
5 changes: 5 additions & 0 deletions PDEInModelica/julia/run.jl
View file
@@ -0,0 +1,5 @@
include("advection.jl")
(X,U) = simulate(0.5)
using Gadfly
plot(x = X, y = U[1,:])
#plot(x = X, y = U[2,:])
74 changes: 74 additions & 0 deletions PDEInModelica/julia/solver.jl
View file
@@ -0,0 +1,74 @@
#numerics setup
nX = 100 #number of nodes


#numerics functions

function updateU_x_LW(X,U,U_x)
for i = 2:(size(X,1)-1)
U_x[i] = (U[i+1] - U[i-1])/(X[i+1] - X[i-1])
end
end

function updateU_t(X,U,U_x,U_t,t)
for i = 1:nU
for j = 2:nX-1
U_t[i,j] = utFun(X[j],U[j],U_x[j],t,i)
end
end
end

function updateU_LW(U,U_t,dt)
function UUpdate(i)
(U[i+1] + U[i-1])/2 + dt*U_t[i]
end
u1new = UUpdate(2)
unew = UUpdate(3)
for i = 4:size(U,2)-1
U[i-2] = u1new
u1new = unew
unew = UUpdate(i)
end
U[size(U,2)-2] = u1new
U[size(U,2)-1] = unew
end


#computation
function simulate(tEnd)


#numerics variables
X = linspace(0.0,L,nX) #coordinate
U = Array(Float64,nU,nX) #state fields u[nVar, nNode]
U_x = Array(Float64,nU,nX) #state fields space derivatives ux[nVar, nNode]
U_t = Array(Float64,nU,nX) #state fields time derivatives ut[nVar, nNode]
t = 0.0 #time
dx = L/(nX-1) #space step
#dt #time step
cfl = 0.5


for i = 1:nU
for j = 1:nX
U[i,j] = initFun(i,X[j])
end
end

while t < tEnd
dt = cfl*maxEigValFun()*dx
for i = 1:nU
U[i,1] = BCFun(i,left,t,X,U)
U[i,nX] = BCFun(i,right,t,X,U)
end
updateU_x_LW(X,U,U_x)
updateU_t(X,U,U_x,U_t,t)
updateU_LW(U,U_t,dt)
t = t + dt
end
(X,U)
end




0 comments on commit 3c23445

Please sign in to comment.