Mercurial > repos > public > sbplib

comparison +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
comparison
equal deleted inserted replaced
833:9f4c45a2d271 834:f1f0bf087e1c
1 %% Function that constructs artificial dissipation operators using undivided differences
2 function D = dissipationOperator(m, order, Hinv, scaling)
3 % TBD: Add or remove D_2 and Dp/Dm?
4 % d2=[1 2 1];
5 % D_2=(diag(ones(m-1,1),-1)-2*diag(ones(m,1),0)+ ...
6 % diag(ones(m-1,1),1));
7 % D_2(1,1:3)=[d2];D_2(m,m-2:m)=[d2];
8 % %Dm
9 % DD_m=(diag(ones(m-1,1),+1)-diag(ones(m,1),0));
10 % DD_m(m,m-1:m)=[-1 1];
11 % %Dp
12 % DD_p=(-diag(ones(m-1,1),-1)+diag(ones(m,1),0));
13 % DD_p(1,1:2)=[-1 1];
14
15 switch order
16 case 1
17 DD_1=(diag(ones(m-1,1),+1)-diag(ones(m,1),0));
18 DD_1(m,m-1:m)=[-1 1];
19 D = DD_2'*DD_2;
20 case 2
21 dd2=0*[1 -2 1];
22 DD_2=(diag(ones(m-1,1),-1)-2*diag(ones(m,1),0)+ ...
23 diag(ones(m-1,1),1));
24 DD_2(1,1:3)=[dd2];DD_2(m,m-2:m)=[dd2];
25 D = DD_2'*DD_2;
26 case 3
27 d3=0*[-1 3 -3 1];
28 DD_3=(-diag(ones(m-2,1),-2)+3*diag(ones(m-1,1),-1)-3*diag(ones(m,1),0)+ ...
29 diag(ones(m-1,1),1));
30 DD_3(1:2,1:4)=[d3;d3];
31 DD_3(m,m-3:m)=[d3];
32 D = DD_3'*DD_3;
33 case 4
34 default_arg('scaling', 1/12);
35 d4=0*[1 -4 6 -4 1];
36 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));
37 DD_4(1:2,1:5)=[d4;d4];DD_4(m-1:m,m-4:m)=[d4;d4];
38 D = DD_4'*DD_4;
39 case 5
40 d5=0*[-1 5 -10 10 -5 1];
41 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));
42 DD_5(1:3,1:6)=[d5;d5;d5];
43 DD_5(m-1:m,m-5:m)=[d5;d5];
44 D = DD_5'*DD_5;
45 case 6
46 default_arg('scaling', 1/60);
47 d6=0*[1 -6 15 -20 15 -6 1];
48 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));
49 DD_6(1:3,1:7)=[d6;d6;d6];DD_6(m-2:m,m-6:m)=[d6;d6;d6];
50 D = DD_6'*DD_6;
51 case 7
52 d7=0*[-1 7 -21 35 -35 21 -7 1];
53 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));
54 DD_7(1:4,1:8)=[d7;d7;d7;d7];
55 DD_7(m-2:m,m-7:m)=[d7;d7;d7];
56 D = DD_7'*DD_7;
57 case 9
58 d9=0*[-1 9 -36 84 -126 126 -84 36 -9 1];
59 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));
60 DD_9(1:5,1:10)=[d9;d9;d9;d9;d9];
61 DD_9(m-3:m,m-9:m)=[d9;d9;d9;d9];
62 D = DD_9'*DD_9;
63 otherwise
64 error('Order not yet supported', order);
65 end
66 scaling
67 D = scaling*Hinv*D;
68 end