Mercurial > repos > public > sbplib

diff +sbp/dissipationOperator.m @ 834:f1f0bf087e1c feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add support for artificial viscosity - Add option to use upwind operators for discretizing the scheme - Add option to use dissipation operators constructed from undivided differences when discretizing with narrow stencil operators
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 13 Sep 2018 14:21:47 +0200
parents
children cda996e64925
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/dissipationOperator.m	Thu Sep 13 14:21:47 2018 +0200
@@ -0,0 +1,68 @@
+%%  Function that constructs artificial dissipation operators using undivided differences
+function D = dissipationOperator(m, order, Hinv, scaling)
+    % TBD: Add or remove D_2 and Dp/Dm?
+    % d2=[1 2 1];
+    % D_2=(diag(ones(m-1,1),-1)-2*diag(ones(m,1),0)+ ...
+    % diag(ones(m-1,1),1));
+    % D_2(1,1:3)=[d2];D_2(m,m-2:m)=[d2];
+    % %Dm
+    % DD_m=(diag(ones(m-1,1),+1)-diag(ones(m,1),0));
+    % DD_m(m,m-1:m)=[-1 1];
+    % %Dp
+    % DD_p=(-diag(ones(m-1,1),-1)+diag(ones(m,1),0));
+    % DD_p(1,1:2)=[-1 1];
+
+    switch order
+        case 1
+            DD_1=(diag(ones(m-1,1),+1)-diag(ones(m,1),0));
+            DD_1(m,m-1:m)=[-1 1];
+            D = DD_2'*DD_2;
+        case 2
+            dd2=0*[1 -2 1];
+            DD_2=(diag(ones(m-1,1),-1)-2*diag(ones(m,1),0)+ ...
+            	  diag(ones(m-1,1),1));
+            DD_2(1,1:3)=[dd2];DD_2(m,m-2:m)=[dd2];
+            D = DD_2'*DD_2;
+        case 3
+            d3=0*[-1 3 -3 1];
+            DD_3=(-diag(ones(m-2,1),-2)+3*diag(ones(m-1,1),-1)-3*diag(ones(m,1),0)+ ...
+                  diag(ones(m-1,1),1));
+            DD_3(1:2,1:4)=[d3;d3];
+            DD_3(m,m-3:m)=[d3];
+            D = DD_3'*DD_3;
+        case 4
+            default_arg('scaling', 1/12);
+            d4=0*[1 -4 6 -4 1];
+            DD_4=(diag(ones(m-2,1),2)-4*diag(ones(m-1,1),1)+6*diag(ones(m,1),0)-4*diag(ones(m-1,1),-1)+diag(ones(m-2,1),-2));
+            DD_4(1:2,1:5)=[d4;d4];DD_4(m-1:m,m-4:m)=[d4;d4];
+            D = DD_4'*DD_4;
+        case 5
+            d5=0*[-1 5 -10 10 -5 1];
+            DD_5=(-diag(ones(m-3,1),-3)+5*diag(ones(m-2,1),-2)-10*diag(ones(m-1,1),-1)+10*diag(ones(m,1),0)-5*diag(ones(m-1,1),1)+diag(ones(m-2,1),2));
+            DD_5(1:3,1:6)=[d5;d5;d5];
+            DD_5(m-1:m,m-5:m)=[d5;d5];
+            D = DD_5'*DD_5; 
+        case 6
+            default_arg('scaling', 1/60);
+            d6=0*[1 -6 15 -20 15 -6 1];
+            DD_6=(diag(ones(m-3,1),3)-6*diag(ones(m-2,1),2)+15*diag(ones(m-1,1),1)-20*diag(ones(m,1),0)+15*diag(ones(m-1,1),-1)-6*diag(ones(m-2,1),-2)+diag(ones(m-3,1),-3));
+            DD_6(1:3,1:7)=[d6;d6;d6];DD_6(m-2:m,m-6:m)=[d6;d6;d6];
+            D = DD_6'*DD_6;
+        case 7
+            d7=0*[-1 7 -21 35 -35 21 -7 1]; 
+            DD_7=(-diag(ones(m-4,1),-4)+7*diag(ones(m-3,1),-3)-21*diag(ones(m-2,1),-2)+35*diag(ones(m-1,1),-1)-35*diag(ones(m,1),0)+21*diag(ones(m-1,1),1)-7*diag(ones(m-2,1),2)+diag(ones(m-3,1),3));
+            DD_7(1:4,1:8)=[d7;d7;d7;d7];
+            DD_7(m-2:m,m-7:m)=[d7;d7;d7];
+            D = DD_7'*DD_7;
+        case 9
+            d9=0*[-1 9 -36 84 -126 126 -84 36 -9 1]; 
+            DD_9=(-diag(ones(m-5,1),-5)+9*diag(ones(m-4,1),-4)-36*diag(ones(m-3,1),-3)+84*diag(ones(m-2,1),-2)-126*diag(ones(m-1,1),-1)+126*diag(ones(m,1),0)-84*diag(ones(m-1,1),1)+36*diag(ones(m-2,1),2)-9*diag(ones(m-3,1),3)+diag(ones(m-4,1),4));
+            DD_9(1:5,1:10)=[d9;d9;d9;d9;d9];
+            DD_9(m-3:m,m-9:m)=[d9;d9;d9;d9];
+            D = DD_9'*DD_9;
+        otherwise
+            error('Order not yet supported', order);
+    end
+    scaling
+    D = scaling*Hinv*D;
+end