PDE extension

experimental string model finished, arterial model started

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

PDEInModelica/doc/modelExamples.lyx

@@ -304,6 +304,19 @@ or
 \end_inset
 
 
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula 
+\begin{multline*}
+\frac{\partial^{2}u}{\partial t^{2}}-c\frac{\partial^{2}v}{\partial t\partial x}=\\
+=\frac{\partial^{2}u}{\partial t^{2}}-c\frac{\partial^{2}v}{\partial x\partial t}=\\
+\frac{\partial^{2}u}{\partial t^{2}}-c\frac{\partial^{2}u}{\partial x^{2}}
+\end{multline*}
+
+\end_inset
+
+
 \end_layout
 
 \begin_layout Paragraph

PDEInModelica/julia/artery.jl

@@ -0,0 +1,76 @@
+#model:
+#  A_t + A_x*U + A*U_x = 0
+#  U_t + (2*alpha-1)*U*U_x + (alpha-1)*U*U/A*A_x + 1/rho*P_x = f/(rho*A)
+#  f = -2*(zeta+2)*mu*C.pi*U
+#  P = P_ext + beta/A_0*(sqrt(A) - sqrt(A_0))
+#  Q = Q_heart
+#  Q = P/R_out
+
+include("incl.jl")
+
+# - model static data
+nU = 2                 #number of state fields u[1,:]
+nV = 2                 #number of algebraic fields
+
+# - model values
+
+L = 2.0                #length of domain
+
+type Parameters
+    alpha
+    zeta
+    rho
+    mu
+    P_ext
+    A_0
+    beta
+    h
+    E
+    CO
+    MAP
+    R_out
+end
+
+p = Parameters(1.1, 0.0, 1000, 4e-3, 0.0, 24e-3, 0.0, 0.002, 6500000.0, 5.6/1000/6, 90*133.322387415, 0.0)
+
+function initializeBoundParameters()
+    p.zeta = (2 - p.alpha)/(p.alpha-1)
+    p.beta = 4.0/3.0*sqrt(pi)*h*E
+end
+
+
+#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 utFun(x,u,ux,t,)
+    [c*ux[2]; ux[1]]
+end
+
+include("solver.jl")

PDEInModelica/julia/incl.jl

@@ -1,7 +1,7 @@
 const left = false
 const right = true
 
-function extrapolate(X,U,nState,side)
+function extrapolate(nState,side,X,U)
     if side == left 
         extrapolate3(X[1],X[2],X[3],X[4],U[nState,2],U[nState,3],U[nState,4])
     elseif side == right

PDEInModelica/julia/plot.p

@@ -1,2 +1 @@
-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,6 +1,6 @@
 #include("advection.jl")
 include("string.jl")
-(X,U) = simulate(0.5)
+(X,U) = simulate(1.0)
 
 function writeData(X,U)
 f = open("result.txt","w")

PDEInModelica/julia/solver.jl

@@ -1,36 +1,34 @@
 #numerics setup
-nX = 400               #number of nodes
+nX = 6400               #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])
+        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
+    for j = 2:nX-1 
+        U_t[:,j] = utFun(X[j],U[:,j],U_x[:,j],t)
     end
 end
 
 function updateU_LF(U,U_t,dt)
     function UUpdate(i)
-        (U[i+1] + U[i-1])/2 + dt*U_t[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
+        U[:,i-2] = u1new
+        u1new = unew[:,:]
         unew = UUpdate(i)
     end
-    U[size(U,2)-2] = u1new
-    U[size(U,2)-1] = unew
+    U[:,size(U,2)-2] = u1new
+    U[:,size(U,2)-1] = unew
 end
 
 
@@ -40,13 +38,14 @@ 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]
+    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
-    cfl = 0.01
+    cfl = 0.2
 
 
     for i = 1:nU
@@ -62,6 +61,8 @@ function simulate(tEnd)
             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)
         t = t + dt

PDEInModelica/julia/strTest.jl

@@ -0,0 +1,44 @@
+include("string.jl")
+tEnd = 1
+(X,U) = simulate(tEnd)
+nX = size(X,1)
+function analyticSol(x,t)
+    v = sqrt(c)
+    if x >= c*t 
+        initFun(1, x - v*t)
+    else
+        l1BCFun(t - x/v)
+    end
+end
+A = Array(Float64,nX)
+for i in 1:nX
+    A[i] = analyticSol(X[i],tEnd)
+end
+
+
+function writeData(X,A,U)
+f = open("result.txt","w")
+(nU,nX) = size(U)
+    for iX = 1:nX
+        write(f,string(X[iX])" ")            
+        write(f,string(A[iX])" ")            
+        for iU = 1:nU
+            write(f,string(U[iU,iX])" ")
+        end
+        write(f,"\n")
+    end
+    close(f)
+end
+
+writeData(X,A,U)
+
+#using Gadfly
+#using Winston
+#plot(x = X, y = U[1,:])
+
+
+D = vec(U[1,:]) - A
+l2err = sqrt((transpose(D)*D)[1])/nX
+println(l2err)
+#plot(x = X, y = D)
+#plot(x = X, y = analyticU)

PDEInModelica/julia/string.jl

@@ -6,7 +6,9 @@ include("incl.jl")
 nU = 2                 #number of state fields u[1,:]
 
 # - model values
-L = 1.0                #length of domain
+L = 2.0                #length of domain
+
+c = 2.0
 
 
 #model functions:
@@ -20,7 +22,7 @@ function initFun(i,x)
 end
 
 function l1BCFun(t)
-    if 0.0 < t < 1.0 sin(2.0*pi*t) else 0.0 end
+    if 0.0 < t < 0.5 sin(2.0*pi*t) else 0.0 end
 end
 
 function BCFun(nState,side,t,X,U)
@@ -31,22 +33,16 @@ function BCFun(nState,side,t,X,U)
             0.0
         end
     else
-        extrapolate(X,U,nState,side)
+        extrapolate(nState,side,X,U)
     end
 end
 
 function maxEigValFun()
-    2
+    c
 end
 
-function utFun(x,u,ux,t,iState)
-    if iState == 1 #u
-        2*ux[2] 
-    elseif iState == 2 #v
-        ux[1]
-    else
-        error("wrong variable index in utFun()")
-    end
+function utFun(x,u,ux,t,)
+    [c*ux[2]; ux[1]]
 end
 
 include("solver.jl")

PDEInModelica/models/arterialPulsWave_discretized.mo

@@ -0,0 +1,51 @@
+model arterialPulsWave
+  constant Integer N=120;
+  constant Real L=4;
+  constant Real dx=L/(N - 1);
+  constant Real x[N]=array(dx*i for i in 1:N);
+  constant Real Pa2mmHgR=133.322387415;
+  parameter Real rho=1000;
+  parameter Real f=0;
+  parameter Real alpha=0;
+  parameter Real h0=0.002;
+  parameter Real E=6500000000.0;
+  parameter Real nu=1/2;
+  parameter Real beta=sqrt(Modelica.Constants.pi)*h0*E/((1 - nu^2)*A0);
+  parameter Real Pext=0;
+  parameter Real A0=Modelica.Constants.pi*0.012^2;
+  parameter Real HR=70/60;
+  parameter Real Tc=1/HR;
+  parameter Real MAP=90*Pa2mmHgR;
+  parameter Real CO=5.6/1000/60;
+  parameter Real SV=CO/HR;
+  parameter Real Qmax=3*Modelica.Constants.pi*SV/(2*Tc);
+  parameter Real Rout=MAP/CO;
+  parameter Real AInit=((MAP - Pext)/beta + sqrt(A0))^2;
+  Real A[N](each start=AInit, each fixed=true);
+  Real Q[N];
+  Real u[N];
+  Real P[N];
+  Real tp;
+//initial conditions:
+initial equation
+  for i in 2:N - 1 loop
+    Q[i]=CO;
+  end for;
+equation
+//border conditions:
+  tp=mod(time, Tc);
+  Q[1]=if tp < Tc/3 then Qmax*sin(3*Modelica.Constants.pi*tp/Tc)^2 else 0;
+  A_x[1]=(A[2] - A[1])/dx;
+  Q_x[1]=(Q[2] - Q[1])/dx;
+
+  Q[N]=0;
+  A_x[N]=(A[N] - A[N - 1])/dx;
+  Q_x[N]=(Q[N] - Q[N - 1])/dx;
+
+//equations
+  pder(A,t) + pder(Q,x) = 0 ??in omega;
+  pder(Q,t) + alpha*(2*Q*pder(Q,x)/A - Q^2*pder(A,x)/A^2) + A/rho*pder(P,x) = f/rho ?? in omega;
+  P = Pext + beta*(sqrt(A) - sqrt(A0));
+//  pder(P,x) = beta/2*pder(A,x)/sqrt(A); - generate automaticaly by deriving the above equation
+  u = Q/A;
+end arterialPulsWave;