PDE tests

started some 2D models git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk/doc@24553 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
OpenModelica · Feb 12, 2015 · f4afd92 · f4afd92
1 parent c7a4740
commit f4afd92
Show file tree

Hide file tree

Showing 7 changed files with 390 additions and 49 deletions.
diff --git a/PDEInModelica/doc/modelExamples.lyx b/PDEInModelica/doc/modelExamples.lyx
@@ -1854,6 +1854,10 @@ key "Horacek"
 
 \end_layout
 
+\begin_layout Subsection
+Signals in axons
+\end_layout
+
 \begin_layout Subsection
 Vibrating membrane (drum) in air
 \begin_inset CommandInset label
@@ -2094,6 +2098,51 @@ lstparams "breaklines=true,caption={Vibrating membrane in air},captionpos=b,fram
 \end_inset
 
 
+\end_layout
+
+\begin_layout Subsection
+Conservation laws in general
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula 
+\[
+\frac{\partial w}{\partial t}+\nabla\cdot F(w)=0
+\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+i.e.
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula 
+\[
+\frac{\partial w}{\partial t}+\left(F(w)\right)=0
+\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Subsection
+Shallow water equations
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula 
+\begin{eqnarray*}
+\frac{\partial\eta}{\partial t}+\nabla\cdot(\eta\cdot\bar{u}) & = & 0\\
+\end{eqnarray*}
+
+\end_inset
+
+
 \end_layout
 
 \begin_layout Subsection
@@ -2316,6 +2365,66 @@ p_{e} & =\varepsilon_{e}\varrho(\gamma-1) & p_{i} & =\varepsilon_{i}\varrho(\gam
 \end_inset
 
 
+\end_layout
+
+\begin_layout Section
+2D examples
+\end_layout
+
+\begin_layout Subsection
+Advection
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula 
+\[
+\frac{\partial\psi}{\partial t}+\nabla\cdot\left(\mathbf{u}\text{\psi}\right)=0
+\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+i.e.
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula 
+\[
+\frac{\partial\psi}{\partial t}+u_{x}\frac{\partial\psi}{\partial x}+u_{y}\frac{\partial\psi}{\partial y}=0
+\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Subsection
+Conservation laws
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula 
+\[
+\frac{\partial w}{\partial t}+\nabla\cdot F(w)=0
+\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula 
+\[
+\frac{\partial w}{\partial t}+\frac{\partial}{\partial x}\left(F_{1}(w)\right)+\frac{\partial}{\partial y}\left(F_{2}(w)\right)=0
+\]
+
+\end_inset
+
+
 \end_layout
 
 \end_body

diff --git a/PDEInModelica/models/advection.mo b/PDEInModelica/models/advection.mo
@@ -1,15 +1,15 @@
 model advection "advection equation"
-  annotation(experiment(GridNodes = 100));
-
   parameter Real L = 1; // length
-
-  parameter PDEDomains.DomainLineSegment1D omega(l = L);
-
-  field Real u(domain = omega, start = 0);
-
+  parameter DomainLineSegment1D omega(l = L);
+  field Real u(domain = omega);
   parameter Real c = 1;
-
+initial equation
+  u = if omega.x<0.25 then cos(2*3.14*omega.x) else 0 indomain omega;
 equation
-  der(u) + c*der(u,x) = 0;        //by default in omega.interior
-  u = sin(2*3.14*time)            in omega.left;
+  der(u) + c*pder(u,x) = 0 indomain omega;
+  u = 1                    indomain omega.left;
+  annotation(experiment(GridNodes = 100));
 end advection;
+
+
+
diff --git a/PDEInModelica/models/advectionDiscretized.mo b/PDEInModelica/models/advectionDiscretized.mo
@@ -1,37 +1,25 @@
 model advectionDiscretized
   // u_t + u_x = 0
-  constant Integer N = 100;
-
   parameter Real L = 1;
-
+  constant Integer N = 100;
   parameter Real dx = L / (N - 1);
   parameter Real[N] x = array(i * dx for i in 0:N - 1);
-
   Real[N] u, u_x;
-
   parameter Real c = 1;
-
 initial equation
-  for i in 2:N - 1 loop
-    //initial conditions:
-    u[i] = 0;
+  for i in 1:N loop
+    u[i] = if x[i]<0.25 then cos(2*3.14*x[i]) else 0;
   end for;
-
 equation
   //unused array elements, eqs. just for balanced system:
-  u_x[1] = 0;
-  u_x[N] = 0;
+  u_x[1] = 0; u_x[N] = 0;
   for i in 2:N - 1 loop
     //discretization of spatial derivative:
-    u_x[i] = (u[i + 1] - u[i - 1]) / dx;
+    u_x[i] = (u[i + 1] - u[i - 1]) / (2*dx);
     // the equation:
     der(u[i]) + c*u_x[i] = 0;
   end for;
-
-  //left BC:
-  u[1] = sin(2 * 3.14 * time);
-
-  //extrapolation in the last node
-  u[N] = 2 * u[N - 1] - u[N - 2];
-
+  u[1] = 1;  //left BC
+  u[N] = 2 * u[N - 1] - u[N - 2]; //extrapolation in the last node
+  annotation(experiment(Interval = 0.002));
 end advectionDiscretized;
diff --git a/PDEInModelica/models/tests/advection2D.mo b/PDEInModelica/models/tests/advection2D.mo
@@ -0,0 +1,67 @@
+package advection2D
+  model advectionBase
+    constant Integer Nx = 40, Ny = 40;
+    parameter Real Lx = 10, Ly = 10;
+    parameter Real dx = Lx / (Nx - 1), dy = Ly / (Ny - 1);
+    parameter Real[Nx] x = array(i * dx for i in 0:Nx - 1);
+    parameter Real[Ny] y = array(i * dy for i in 0:Ny - 1);
+    Real[Nx, Ny] u;
+    Real[Nx, Ny] u_x;
+    Real[Nx, Ny] u_y;
+    parameter Real ax = 1, ay = 1;
+  equation
+    //boundary:
+    ///dummy:
+    u_x[1, :] = zeros(Ny);
+    u_x[Nx, :] = zeros(Ny);
+    u_x[2:Nx - 1, 1] = zeros(Nx - 2);
+    u_x[2:Nx - 1, Ny] = zeros(Nx - 2);
+    u_y[1, :] = zeros(Ny);
+    u_y[Nx, :] = zeros(Ny);
+    u_y[2:Nx - 1, 1] = zeros(Nx - 2);
+    u_y[2:Nx - 1, Ny] = zeros(Nx - 2);
+    //interior:
+    for i in 2:Nx - 1 loop
+      for j in 2:Ny - 1 loop
+        //differences
+        u_x[i, j] = (u[i + 1, j] - u[i - 1, j]) / (2 * dx);
+        u_y[i, j] = (u[i, j + 1] - u[i, j - 1]) / (2 * dy);
+        //equations
+        der(u[i, j]) + ax * u_x[i, j] + ay * u_y[i, j] = 0;
+      end for;
+    end for;
+  end advectionBase;
+
+  function icfun
+    input Real x, y;
+    output Real u;
+  protected
+  algorithm
+    u := 0;
+  end icfun;
+
+  model advectionZero
+    extends advectionBase;
+  initial equation
+    u[2:Nx - 1, 2:Ny - 1] = zeros(Nx - 2, Ny - 2);
+  equation
+    u[1, :] = zeros(Ny);
+    u[Nx, :] = zeros(Ny);
+    u[2:Nx - 1, 1] = zeros(Nx - 2);
+    u[2:Nx - 1, Ny] = zeros(Nx - 2);
+  end advectionZero;
+
+  model advectionCosPeak
+    extends advectionBase;
+  initial equation
+    for i in 2:Nx - 1, j in 2:Ny - 1 loop
+      u[i, j] = if sqrt((x[i] - 2) ^ 2 + (y[j] - 2) ^ 2) < 1 then cos(Modelica.Constants.pi / 2 * sqrt((x[i] - 2) ^ 2 + (y[j] - 2) ^ 2)) else 0;
+    end for;
+  equation
+    u[1, :] = zeros(Ny);
+    u[Nx, :] = zeros(Ny);
+    u[2:Nx - 1, 1] = zeros(Nx - 2);
+    u[2:Nx - 1, Ny] = zeros(Nx - 2);
+    annotation(experiment(StartTime = 0, StopTime = 0.2, Tolerance = 1, Interval = 0.05));
+  end advectionCosPeak;
+end advection2D;
diff --git a/...Modelica/models/tests/conservationLaws.mo → ...delica/models/tests/conservationLaws1D.mo b/...Modelica/models/tests/conservationLaws.mo → ...delica/models/tests/conservationLaws1D.mo
@@ -1,4 +1,4 @@
-package conservationLaws
+package conservationLaws1D
   model conservationLaw
     constant Integer M;
     constant Integer N;
@@ -15,7 +15,7 @@ package conservationLaws
   end pokusy;
 
   model advection
-    extends conservationLawLF;
+    extends conservationLawCentral;
     Real u[N];
     parameter Real a = 1;
   initial equation
@@ -93,34 +93,34 @@ package conservationLaws
   end eulerEq;
 
   model Riemann1
-    extends eulerEq(rho_l = 1, u_l = 0.75, p_l = 1, rho_r = 0.125, u_r = 0, p_r = 0.1, x_0 = 0.3, N = 1000, dt = 2e-4, alpha = 0.01);
+    extends eulerEq(rho_l = 1, u_l = 0.75, p_l = 1, rho_r = 0.125, u_r = 0, p_r = 0.1, x_0 = 0.3, N = 1000, dt = 0.0002, alpha = 0.4);
     //euler, N = 200 (dx = 0.005, dt = 0.0001, alpha = 0.01 - celkem maká, hodně difuse, trochu osciluje
     //euler, N = 1000, dt = 2e-5, alpha = 0.01 - dobrý, malinko osciluje na čele první vlny
     //radau1, N = 100, dt = 0.002, alpha = 0.01 - vypadá nejlíp
-    annotation(experiment(StartTime = 0, StopTime = 0.2, Tolerance = 1e-6, Interval = 2e-4));
+    annotation(experiment(StartTime = 0, StopTime = 0.2, Tolerance = 1, Interval = 0.0002));
   end Riemann1;
 
   model advectionCos
-    extends advection(M = 1, N = 100, dt = 0.01, alpha = 1);
+    extends advection(M = 1, N = 1600, dt = 0.0003125);
   initial equation
     u = array(if x[i] < 0.25 then cos(Modelica.Constants.pi / 2 * 4 * x[i]) else 0 for i in 1:N);
   equation
     //BC:
     u[1] = 1;
     u[N] = 0;
-    annotation(experiment(StartTime = 0, StopTime = 0.4, Tolerance = 1e-06, Interval = 0.01));
+    annotation(experiment(StartTime = 0, StopTime = 0.4, Tolerance = 1, Interval = 0.0003125));
   end advectionCos;
 
   model advectionStep
-    extends advection(M = 1, N = 100, dt = 0.009, alpha = 1);
+    extends advection(M = 1, N = 400, dt = 0.0005);
     parameter Real ul = 1, ur = 0;
   initial equation
     u = array(if x[i] < 0.2 then ul else ur for i in 1:N);
   equation
     //BC:
     u[1] = ul;
     u[N] = ur;
-    annotation(experiment(StartTime = 0, StopTime = 0.4, Tolerance = 1e-06, Interval = 0.009));
+    annotation(experiment(StartTime = 0, StopTime = 0.4, Tolerance = 1, Interval = 0.0005));
   end advectionStep;
   annotation(Icon(coordinateSystem(extent = {{-100, -100}, {100, 100}}, preserveAspectRatio = true, initialScale = 0.1, grid = {2, 2})), Diagram(coordinateSystem(extent = {{-100, -100}, {100, 100}}, preserveAspectRatio = true, initialScale = 0.1, grid = {2, 2})));
-end conservationLaws;
+end conservationLaws1D;