PDE extension experiment

Plotting scripts

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

PDEInModelica/models/tests/conservationLaws.mo

@@ -15,26 +15,21 @@ package conservationLaws
   end pokusy;
 
   model advection
-    extends conservationLawLF(M = 1, N = 100, dt = 0.005);
+    extends conservationLawLF;
     Real u[N];
-    parameter Real a = 1, ul = 1, ur = 0;
+    parameter Real a = 1;
   initial equation
-    u = array(if i < N / 2 then ul else ur for i in 1:N);
     //cos BC
     //array(if x[i] < 0.25 then cos(Modelica.Constants.pi / 2 * 4 * x[i]) else 0 for i in 1:N);
     //step BC
-    //array(if i < N / 2 then ul else ur for i in 1:N);
+    //
   equation
-    //BCs
-    u[1] = ul;
-    u[N] = ur;
     //dummy BCs
     F_x[1, 1] = 0;
     F_x[1, N] = 0;
     //equations
     W[1, :] = u;
     F[1, :] = a * u;
-    annotation(experiment(StartTime = 0, StopTime = 0.2, Tolerance = 1e-06, Interval = 0.005));
   end advection;
 
   model conservationLawCentral
@@ -51,20 +46,22 @@ package conservationLaws
 
   model conservationLawLF
     parameter Real dt;
+    parameter Real alpha = 1;
     extends conservationLaw;
   equation
     for i in 2:N - 1 loop
       //discretization of spatial derivative:
       F_x[:, i] = (F[:, i + 1] - F[:, i - 1]) / (2 * dx);
       // the equation:
-      der(W[:, i]) - (W[:, i + 1] - 2 * W[:, i] + W[:, i - 1]) / (2 * dt) + F_x[:, i] = zeros(M);
+      der(W[:, i]) - alpha * (W[:, i + 1] - 2 * W[:, i] + W[:, i - 1]) / (2 * dt) + F_x[:, i] = zeros(M);
       /*- (W[:, i + 1] - 2 * W[:, i] + W[:, i - 1]) / (2 * dt)*/
     end for;
   end conservationLawLF;
 
   model eulerEq
-    extends conservationLawCentral(M = 3);
+    extends conservationLawLF(M = 3);
     Real[N] rho, u, p, E;
+    //	parameter Real u_s;
     parameter Real gamma = 1.4, x_0, T;
     parameter Real rho_l, u_l, p_l;
     parameter Real rho_r, u_r, p_r;
@@ -92,11 +89,38 @@ package conservationLaws
     F[3, :] = u .* (E + p);
     //state equation
     E = p ./ (gamma - 1.0) + rho .* u .^ 2 / 2.0;
+    //  u_s = sqrt(gamma * p ./ rho);
   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 = 0.0005);
-    annotation(experiment(StartTime = 0, StopTime = 0.2, Tolerance = 1e-06, Interval = 0.0005));
+    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);
+    //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));
   end Riemann1;
+
+  model advectionCos
+    extends advection(M = 1, N = 100, dt = 0.01, alpha = 1);
+  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));
+  end advectionCos;
+
+  model advectionStep
+    extends advection(M = 1, N = 100, dt = 0.009, alpha = 1);
+    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));
+  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;

PDEInModelica/models/tests/eulerTests.mo

@@ -1,11 +1,12 @@
 package eulerTests
   model Euler
-    constant Integer N = 100;
+    constant Integer N = 800;
     parameter Real L = 1;
     parameter Real dx = L / (N - 1);
+    parameter Real dt;
     parameter Real[N] x = array(i * dx for i in 0:N - 1);
-    Real[N] rho, v, p, eps;
-    Real[N] rho_x, v_x, p_x;
+    Real[N] rho, v, p, E;
+    Real[N] rho_x, v_x, p_x, E_x;
     parameter Real gamma = 1.4, x_0, T;
     parameter Real rho_l, v_l, p_l;
     parameter Real rho_r, v_r, p_r;
@@ -21,18 +22,22 @@ package eulerTests
     v_x[N] = 0;
     p_x[1] = 0;
     p_x[N] = 0;
+    E_x[1] = 0;
+    E_x[N] = 0;
     for i in 2:N - 1 loop
       //discretization of spatial derivative:
       rho_x[i] = (rho[i + 1] - rho[i - 1]) / (2 * dx);
       v_x[i] = (v[i + 1] - v[i - 1]) / (2 * dx);
       p_x[i] = (p[i + 1] - p[i - 1]) / (2 * dx);
+      E_x[i] = (E[i + 1] - E[i - 1]) / (2 * dx);
     end for;
     // the equation:
     der(rho) = (-rho_x .* v) - rho .* v_x;
     der(rho .* v) = (-rho_x .* v .^ 2) - rho .* v .* v_x - p_x;
-    rho .* der(eps) = -p .* v_x;
-    p = eps .* rho .* (gamma - 1);
-    annotation(experiment(StartTime = 0, StopTime = 0.2, Tolerance = 1e-08, Interval = 2e-05));
+    der(E) = (-v_x .* (E + p)) - v .* (E_x + p_x);
+    //-v_x .* (E + p) - v .* (E_x + p_x);
+    E = p / (gamma - 1) + rho .* v .^ 2 / 2.0;
+    annotation(experiment(StartTime = 0, StopTime = 0.2, Tolerance = 1e-08, Interval = 0.001));
   end Euler;
 
   model EulerViscosity
@@ -79,8 +84,8 @@ package eulerTests
   end EulerViscosity;
 
   model Riemann1
-    Euler riemann1(rho_l = 1, v_l = 0.75, p_l = 1, rho_r = 0.125, v_r = 0, p_r = 0.1, x_0 = 0.3);
-    annotation(experiment(StartTime = 0, StopTime = 0.2, Tolerance = 1, Interval = 0.0002));
+    extends Euler(rho_l = 1, v_l = 0.75, p_l = 1, rho_r = 0.125, v_r = 0, p_r = 0.1, x_0 = 0.3);
+    annotation(experiment(StartTime = 0, StopTime = 0.2, Tolerance = 1, Interval = 0.001));
   end Riemann1;
 
   model Riemann1V

PDEInModelica/models/tests/plotMat.m

@@ -1,19 +1,53 @@
-%fileName = "eulerTests.Riemann1V_res.mat";
-fileName = "conservationLaws.advection_res.mat";
+%fileName = "eulerTests.Riemann1V";
+%fileName = "conservationLaws.advectionCos";
+fileName = "conservationLaws.advectionStep";
 varName = "u[";
-load(fileName);
+xName = "x[";
+exactSol = str2func(strrep(fileName,'.','_'));
+load([fileName "_res.mat"]);
 nameT = transpose(name);
-indexes = [];
+indexesV = [];
+indexesX = [];
 for i = 1:size(name,2)
   if isequal(varName,nameT(i,1:size(varName,2)))
-	  indexes = [indexes i];
+	  indexesV = [indexesV i];
   end;
+  if isequal(xName,nameT(i,1:size(xName,2)))
+	  indexesX = [indexesX i];
+  end;
+end;
+if size(indexesV) != size(indexesX) 
+  error("different sizes of X and variable")
+end;
+size(dataInfo);
+ar = dataInfo(1,indexesX) != 1;
+if ar
+  error("some X isnt in data_1")
+end;
+ar = dataInfo(1,indexesV) != 2;
+if ar
+  error("some var isnt in data_2")
 end;
+indexesX1 = dataInfo(2,indexesX);
+indexesV1 = dataInfo(2,indexesV);
+size(data_1)
+X = data_1(indexesX1,1);
 nSteps = size(data_2,2)
-nPlots = 10;
+nPlots = 5;
+h = figure();
+
+
+
 for i = 1:nPlots
- nFrame = round(1 + (nSteps-1)/(nPlots-1)*(i-1) ) 
- var = data_2(indexes,nFrame);
- plot(var);
- pause(1);
+ nFrame = round(1 + (nSteps-1)/(nPlots-1)*(i-1) ) ;
+ time = data_2(1,nFrame);
+ var = data_2(indexesV,nFrame);
+ varE = feval(exactSol,time,X);
+ plot(X,var,X,varE);
+ pause(0.1);
 end;
+saveas(h,[fileName ".png"],"png");
+m = max(varE);
+err = sum(abs((varE - var)./max(m/100,varE).*(varE>m/100)))
+
+

PDEInModelica/models/tests/plotMat2.m

@@ -1,11 +1,12 @@
-fileName = "conservationLaws.Riemann1_res.mat";
-%fileName = "conservationLaws.advection_res.mat";
-varName = "rho";
+%fileName = "eulerTests.Riemann1_res.mat";
+%fileName = "conservationLaws.Riemann1_res.mat";
+fileName = "conservationLaws.advection_res.mat";
+varName = "u";
 [nameT, data_2, N] = openMat(fileName);
 nPlots = 10;
 for i = 1:nPlots
  nFrame = round(1 + (N-1)/(nPlots-1)*(i-1) ) ;
  getTimeMat(nFrame, nameT, data_2) 
  plot(getVarMat(varName, nFrame, nameT, data_2));
- pause(1);
+ pause(0.2);
 end;