PDE extension experimental solver

working on algebraic variables

git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk/doc@22399 f25d12d1-65f4-0310-ae8a-bbce733d8d8e

Author: Jan Silar

PDEInModelica/julia/plot.p

@@ -1 +1,2 @@
-plot "result.txt" using 1:2, "result.txt" using 1:3
+plot "result.txt" using 1:2
+#, "result.txt" using 1:3

PDEInModelica/julia/run.jl

@@ -1,5 +1,6 @@
 #include("advection.jl")
 include("string.jl")
+#include("stringAlg.jl")
 (X,U) = simulate(1.0)
 
 function writeData(X,U)

PDEInModelica/julia/solver.jl

@@ -1,28 +1,35 @@
 #numerics setup
-nX = 6400               #number of nodes
+nX = 400               #number of nodes
 
 
 #numerics functions
 
-function updateU_x_LF(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])
+function diff_x_LF(X,Y,Y_x)
+    for i = 2:(size(X,1)-1) #over inner grid nodes
+        Y_x[:,i] = (Y[:,i+1] - Y[:,i-1])/(X[i+1] - X[i-1])
     end
 end
 
-function updateU_t(X,U,U_x,U_t,t)
-    for j = 2:nX-1 
-        U_t[:,j] = utFun(X[j],U[:,j],U_x[:,j],t)
+function updateV(X,U,U_x,t,V)
+    for i = 1:size(X,1) #over all grid nodes
+#FIX: V[:,1] and V[:,nX] are evaluated using unassigned U_x
+        V[:,i] = vFun(X[i],U[:,i],U_x[:,i],t)
     end
 end
 
-function updateU_LF(U,U_t,dt)
+function updateU_t(X,U,U_x,V,V_x,t,U_t)
+    for i = 2:nX-1 #over inner grid nodes
+        U_t[:,i] = utFun(X[i],U[:,i],U_x[:,i],V[:,i],V_x[:,i],t)
+    end
+end
+
+function updateU_LF(U_t,dt,U)
     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 
+    for i = 4:size(U,2)-1 #over rest of inner grid nodes
         U[:,i-2] = u1new
         u1new = unew[:,:]
         unew = UUpdate(i)
@@ -35,16 +42,17 @@ 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]
-    V = Array(Float64,nV,nX) #algevraic fields  V[nVar, nNode]
-    t = 0.0                #time
-    dx = L/(nX-1)          #space step
-    #dt                     #time step
+    #model 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 derivative Ux[nVar, nNode]
+    U_t = Array(Float64,nU,nX) #state fields time derivative Ut[nVar, nNode]
+    V = Array(Float64,nV,nX)   #algevraic fields  V[nVar, nNode]
+    V_x = Array(Float64,nV,nX) #algevraic fields space derivative V[nVar, nNode]
+    t = 0.0                    #time
+    #numerics variables:
+    dx = L/(nX-1)              #space step
+    #dt                         #time step
     cfl = 0.2
 
 
@@ -60,11 +68,11 @@ function simulate(tEnd)
             U[i,1]  = BCFun(i,left,t,X,U)
             U[i,nX] = BCFun(i,right,t,X,U)
         end
-        updateU_x_LF(X,U,U_x)
-#TODO        updateV
-#TODO        updateV_x
-        updateU_t(X,U,U_x,U_t,t)
-        updateU_LF(U,U_t,dt)
+        diff_x_LF(X,U,U_x)
+        updateV(X,U,U_x,t,V)
+        diff_x_LF(X,V,V_x)
+        updateU_t(X,U,U_x,V,V_x,t,U_t)
+        updateU_LF(U_t,dt,U)
         t = t + dt
     end    
     (X,U)

PDEInModelica/julia/stringAlg.jl

@@ -0,0 +1,54 @@
+#model pder(u,t) + 1*pder(u,x) = 0
+
+include("incl.jl")
+
+# - model static data
+nU = 2                 #number of state fields u[nu,nx]      u,v
+nV = 1                 #number of algebraic fields v[nv,nx]  w
+
+# - model values
+L = 2.0                #length of domain
+
+c = 2.0
+
+
+#model functions:
+
+function initFun(i,x)
+    if i == 1 || i == 2 
+        0.0 
+    else 
+        error("wrong state number in advection") 
+    end 
+end
+
+function l1BCFun(t)
+    if 0.0 < t < 0.5 sin(2.0*pi*t) else 0.0 end
+end
+
+function BCFun(nState,side,t,X,U)
+    if nState == 1
+        if side == left 
+            l1BCFun(t)
+        elseif side == right
+            0.0
+        end
+    else
+        extrapolate(nState,side,X,U)
+    end
+end
+
+function maxEigValFun()
+    c
+end
+
+function vFun(x,u,u_x,t)
+    [c*u[1]]          #w = c*v;
+end
+
+
+function utFun(x,u,ux,v,vx,t)
+    [vx[1]; ux[1]]    #pder(u,time) - pder(w,x) = 0 ;  pder(v,time) - pder(u,x) = 0
+end
+
+include("solver.jl")