PDE extension

Working on artery model git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk/doc@22910 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
OpenModelica · Oct 24, 2014 · e7c59cf · e7c59cf
1 parent 9b26ab2
commit e7c59cf
Show file tree

Hide file tree

Showing 5 changed files with 97 additions and 25 deletions.
diff --git a/PDEInModelica/doc/questionsAndNotes.lyx b/PDEInModelica/doc/questionsAndNotes.lyx
@@ -1117,6 +1117,66 @@ Build whole solver in some PDE framework, perhaps Overture (http://www.overturef
 amework.org/)
 \end_layout
 
+\begin_layout Itemize
+
+\lang british
+libraries to solve PDEs: 
+\end_layout
+
+\begin_deeper
+\begin_layout Itemize
+
+\series bold
+\lang british
+clawpack
+\series default
+, http://www.amath.washington.edu/~claw/
+\end_layout
+
+\begin_layout Itemize
+
+\lang british
+pydstool, http://www.cam.cornell.edu/~rclewley/cgi-bin/moin.cgi/
+\end_layout
+
+\begin_layout Itemize
+
+\lang british
+FemHub, http://femhub.org/welcome.html
+\end_layout
+
+\begin_layout Itemize
+
+\lang british
+Hermes
+\end_layout
+
+\begin_layout Itemize
+
+\series bold
+\lang british
+dealII
+\end_layout
+
+\begin_layout Itemize
+
+\lang british
+overture, http://www.overtureframework.org/
+\end_layout
+
+\begin_layout Itemize
+
+\lang british
+hiflow3, http://www.hiflow3.org/
+\end_layout
+
+\begin_layout Itemize
+
+\lang british
+flex, http://www.pdesolutions.com/
+\end_layout
+
+\end_deeper
 \begin_layout Section
 TODO
 \end_layout

diff --git a/PDEInModelica/julia/artery.jl b/PDEInModelica/julia/artery.jl
@@ -21,7 +21,7 @@ nV = 2                 #number of algebraic fields
 
 L = 2.0                #length of domain
 
-type Parameters
+type P
     alpha
     zeta
     rho
@@ -36,11 +36,12 @@ type Parameters
     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)
+p = P(1.1, 0.0, 1000, 4e-3, 0.0, 24e-3, 0.0, 0.002, 6500000.0, 5.6/1000/60, 90*133.322387415, 0.0)
+#             alpha zeta rho   mu   P_ext A_0    beta h      E          CO           MAP              R_out
 
 function initializeBoundParameters(p)
     p.zeta = (2 - p.alpha)/(p.alpha-1)
-    p.beta = 4.0/3.0*sqrt(pi)*h*E
+    p.beta = 4.0/3.0*sqrt(pi)*p.h*p.E
 end
 
 
@@ -63,36 +64,46 @@ function Q_heart(p,t)
     T1 = mod(t,Tc)
     SV = p.CO/HR
     Q_max = SV / (10 * TcP)
-
-    if T1 < TcP T1/TcP*Q_max
-    elseif T1 < 9*TcP Q_max
-    elseif T1 < 10*TcP (10*TcP-t)/TcP*Q_max
-    else 0.0
+    if T1 < TcP 
+        T1/TcP*Q_max
+    elseif T1 < 9*TcP 
+        Q_max
+    elseif T1 < 10*TcP 
+        (10*TcP-T1)/TcP*Q_max
+    else 
+        0.0
     end
 end
-else 0.0 end
-end
 
-function BCFun(t,X,U)
+function BCFun(p,t,X,U)
     Ql = Q_heart(p,t)
     Al = extrapolate(2,left,X,U)
     Ul = Ql/Al
     Ar = extrapolate(2,right,X,U)
-    Ur = CO/Ar
+    Ur = p.CO/Ar
     ([Ul;Al],[Ur,Ar])
 
 end
 
-function maxEigValFun()
-    c
+function maxEigValFun(p,U,V)
+    sqrt(p.E*p.h/(2*p.rho*sqrt(A/pi))) #TODO: dopsat, upravit signaturu maxEigValFun v solveru a ostatních modelech
 end
 
-function vFun(p,x,u,U_x,t)
-#TODO: implement
+function vFun(p,x,u,ux,t)
+    #P = P_ext + beta/A_0*(sqrt(A) - sqrt(A_0)) :
+    P = p.P_ext + p.beta/p.A_0*(sqrt(u[2]) - sqrt(p.A_0))
+    #f = -2*(zeta+2)*mu*C.pi*U :
+    f = -2*(p.zeta+2)*p.mu*pi*u[1]
+    [P, f]
 end
 
 function utFun(p,x,u,ux,v,vx,t,)
-    [c*ux[2]; ux[1]]
+    #pder(U,time) = f/(rho*A) - (2*alpha-1)*U*pder(U,x) - (alpha-1)*U*U/A*pder(A,x) - 1/rho*pder(P,x)
+    U_t = v[2]/(p.rho*u[2]) - (2*p.alpha-1)*u[1]*ux[1] - (p.alpha-1)*u[1]*u[1]/u[2]ux[2] - 1/p.rho*vx[1]
+    #pder(A,time) = - pder(A*U,x)
+    #                  A_x*U + A*U_x
+    A_t = - ux[2]*u[1] - ux[1]*u[2]
+    [A_t, A_t]
 end
 
 include("solver.jl")
diff --git a/PDEInModelica/julia/plot.p b/PDEInModelica/julia/plot.p
@@ -1,2 +1,2 @@
-plot "result.txt" using 1:2
-#, "result.txt" using 1:3
+#plot "result.txt" using 1:2
+plot "result.txt" using 1:3
diff --git a/PDEInModelica/julia/run.jl b/PDEInModelica/julia/run.jl
@@ -1,7 +1,8 @@
 #include("advection.jl")
 #include("string.jl")
-include("stringAlg.jl")
-(X,U) = simulate(1.0)
+#include("stringAlg.jl")
+include("artery.jl")
+(X,U) = simulate(5.0)
 
 function writeData(X,U)
 f = open("result.txt","w")

diff --git a/PDEInModelica/models/arterialPulsWave.mo b/PDEInModelica/models/arterialPulsWave.mo
@@ -22,11 +22,11 @@ model arterialPulsWave
 initial equation
   U*A = CO                                                      in omega;
 equation
-  pder(A,time) + pder(A*U,x) = 0                                in omega;
+  pder(A,time) + pder(A*U,x) = 0                                in omega;//4 -> A_t (utFun)
   pder(U,time) + (2*alpha-1)*U*pder(U,x) + (alpha-1)*U*U/A*pder(A,x)
-                                  + 1/rho*pder(P,x) = f/(rho*A) in omega;
-  f = -2*(zeta+2)*mu*C.pi*U                                     in omega;
-  P = P_ext + beta/A_0*(sqrt(A) - sqrt(A_0))                    in omega;
+                                  + 1/rho*pder(P,x) = f/(rho*A) in omega;//3 -> U_t (utFun)
+  f = -2*(zeta+2)*mu*C.pi*U                                     in omega;//2 -> f (vFun)
+  P = P_ext + beta/A_0*(sqrt(A) - sqrt(A_0))                    in omega;//1 -> P (vFun)
   A*U = Q_heart;                                          in omega.left;
   A*U = CO;                                               in omega.right;
 end arterialPulsWave;