Mercurial > repos > public > sbplib

changeset 845:1e057b0f2fed feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add RK6 with residual viscosity update and reduce computational effort of spatial scheme - Add RK6 with residual updates - Change the D2 operator for upwind schemes to one less computationally expensive.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 19 Sep 2018 16:32:05 +0200
parents 9e4e0576ca0f
children c6fcee3fcf1b
files +scheme/Burgers1D.m +time/+rk4/rungekutta_6.m +time/+rk4/rungekutta_6RV.m +time/Rungekutta4RV.m
diffstat 4 files changed, 46 insertions(+), 12 deletions(-) [+]
line wrap: on
line diff
diff -r 9e4e0576ca0f -r 1e057b0f2fed +scheme/Burgers1D.m
--- a/+scheme/Burgers1D.m	Tue Sep 18 10:32:00 2018 +0200
+++ b/+scheme/Burgers1D.m	Wed Sep 19 16:32:05 2018 +0200
@@ -37,9 +37,9 @@
                 case 'upwind'
                     ops = sbp.D1Upwind(m, lim, order);
                     D1 = (ops.Dp + ops.Dm)/2;
-                    %D2eps = @(eps) ops.Dp*diag(eps)*ops.Dm;
-                    %D2eps = @(eps) ops.Dm*diag(eps)*ops.Dp;
-                    D2 = @(eps) (ops.Dp*diag(eps)*ops.Dm + ops.Dm*diag(eps)*ops.Dp)/2;
+                    D2 = @(eps) ops.Dp*spdiag(eps)*ops.Dm;
+                    %D2 = @(eps) ops.Dm*diag(eps)*ops.Dp;
+                    %D2 = @(eps) (ops.Dp*diag(eps)*ops.Dm + ops.Dm*diag(eps)*ops.Dp)/2;
                     if (strcmp(dissipation,'on'))
                         DissipationOp = (ops.Dp-ops.Dm)/2;
                     end
diff -r 9e4e0576ca0f -r 1e057b0f2fed +time/+rk4/rungekutta_6.m
--- a/+time/+rk4/rungekutta_6.m	Tue Sep 18 10:32:00 2018 +0200
+++ b/+time/+rk4/rungekutta_6.m	Wed Sep 19 16:32:05 2018 +0200
@@ -1,7 +1,7 @@
-% Takes one time step of size k using the rungekutta method
+% Takes one time step of size dt using the rungekutta method
 % starting from v_0 and where the function F(v,t) gives the
 % time derivatives.
-function v = rungekutta_6(v, t , k, F)
+function v = rungekutta_6(v, t , dt, F)
     s = 7
     k = zeros(length(v),s)
     a = zeros(7,6);
@@ -15,17 +15,17 @@
         229/1200 - 29/6000*sqrt(5),  119/240 - 187/1200*sqrt(5), -14/75 + 34/375*sqrt(5), -3/100*sqrt(5),        0,                 0;
         71/2400 - 587/12000*sqrt(5), 187/480 - 391/2400*sqrt(5), -38/75 + 26/375*sqrt(5), 27/80 - 3/400*sqrt(5), (1+sqrt(5))/4,     0;
         -49/480 + 43/160*sqrt(5),    -425/96 + 51/32*sqrt(5),    52/15 - 4/5*sqrt(5),     -27/16 + 3/16*sqrt(5), 5/4 - 3/4*sqrt(5), 5/2 - 1/2*sqrt(5);
-    ]
+    ];
 
     for i = 1:s
-        u = v
-        for j = 1: i-1
-            u = u + h*a(i,j) * k(:,j)
+        u = v;
+        for j = 1:i-1
+            u = u + dt*a(i,j)*k(:,j);
         end
-        k(:,i) = F(t+c(i)*k,u)
+        k(:,i) = F(t+c(i)*dt,u);
     end
 
     for i = 1:s
-        v = v + k*b(i)*k(:,i)
+        v = v + dt*b(i)*k(:,i);
     end
 end
diff -r 9e4e0576ca0f -r 1e057b0f2fed +time/+rk4/rungekutta_6RV.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk4/rungekutta_6RV.m	Wed Sep 19 16:32:05 2018 +0200
@@ -0,0 +1,34 @@
+% Takes one time step of size dt using the rungekutta method
+% starting from v_0 and where the function F(v,t) gives the
+% time derivatives.
+function v = rungekutta_6RV(v, t , dt, F, RV)
+    s = 7;
+    k = zeros(length(v),s);
+    a = zeros(7,6);
+    c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1];
+    b = [1/12, 0, 0, 0, 5/12, 5/12, 1/12];
+    a = [
+        0,                           0,                          0,                       0,                     0,                 0;
+        4/7,                         0,                          0,                       0,                     0,                 0;
+        115/112,                     -5/16,                      0,                       0,                     0,                 0;
+        589/630,                     5/18,                       -16/45,                  0,                     0,                 0;
+        229/1200 - 29/6000*sqrt(5),  119/240 - 187/1200*sqrt(5), -14/75 + 34/375*sqrt(5), -3/100*sqrt(5),        0,                 0;
+        71/2400 - 587/12000*sqrt(5), 187/480 - 391/2400*sqrt(5), -38/75 + 26/375*sqrt(5), 27/80 - 3/400*sqrt(5), (1+sqrt(5))/4,     0;
+        -49/480 + 43/160*sqrt(5),    -425/96 + 51/32*sqrt(5),    52/15 - 4/5*sqrt(5),     -27/16 + 3/16*sqrt(5), 5/4 - 3/4*sqrt(5), 5/2 - 1/2*sqrt(5);
+    ];
+
+    for i = 1:s
+        u = v;
+        for j = 1:i-1
+            u = u + dt*a(i,j)*k(:,j);
+        end
+        if (i > 1)
+            RV.update(u,v,c(i)*dt);
+        end
+        k(:,i) = F(u,t+c(i)*dt,RV.getViscosity());
+    end
+
+    for i = 1:s
+        v = v + dt*b(i)*k(:,i);
+    end
+end
diff -r 9e4e0576ca0f -r 1e057b0f2fed +time/Rungekutta4RV.m
--- a/+time/Rungekutta4RV.m	Tue Sep 18 10:32:00 2018 +0200
+++ b/+time/Rungekutta4RV.m	Wed Sep 19 16:32:05 2018 +0200
@@ -34,7 +34,7 @@
         end
 
         function obj = step(obj)
-            obj.v = time.rk4.rungekutta_4RV(obj.v, obj.t, obj.k, obj.F, obj.RV);
+            obj.v = time.rk4.rungekutta_6RV(obj.v, obj.t, obj.k, obj.F, obj.RV);
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
         end