Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst3d.m @ 1231:e1f9bedd64a9 feature/volcano

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add opSet to Hypsyst3d
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 18 Nov 2019 17:17:16 -0800
parents ef08adea56c4
children 10881b234f77
comparison
equal deleted inserted replaced
1228:ef08adea56c4 1231:e1f9bedd64a9
11 Aevaluated,Bevaluated,Cevaluated, Eevaluated 11 Aevaluated,Bevaluated,Cevaluated, Eevaluated
12 12
13 H % Discrete norm 13 H % Discrete norm
14 Hx, Hy, Hz % Norms in the x, y and z directions 14 Hx, Hy, Hz % Norms in the x, y and z directions
15 Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. 15 Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
16 I_x,I_y, I_z, I_N
17 e_w, e_e, e_s, e_n, e_b, e_t 16 e_w, e_e, e_s, e_n, e_b, e_t
18 params % Parameters for the coefficient matrices 17 params % Parameters for the coefficient matrices
19 end 18 end
20 19
21 20
22 methods 21 methods
23 % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Cu_z-Eu 22 % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Cu_z-Eu
24 function obj = Hypsyst3d(cartesianGrid, order, A, B, C, E, params, operator) 23 function obj = Hypsyst3d(cartesianGrid, order, A, B, C, E, params, opSet)
25 assertType(cartesianGrid, 'grid.Cartesian'); 24 assertType(cartesianGrid, 'grid.Cartesian');
26 default_arg('E', []) 25 default_arg('E', [])
26 default_arg('opSet', @sbp.D2Standard)
27 obj.grid = cartesianGrid; 27 obj.grid = cartesianGrid;
28 xlim = obj.grid.lim{1};
29 ylim = obj.grid.lim{2};
30 zlim = obj.grid.lim{3};
31 28
32 obj.A = A; 29 obj.A = A;
33 obj.B = B; 30 obj.B = B;
34 obj.C = C; 31 obj.C = C;
35 obj.E = E; 32 obj.E = E;
36 m_x = obj.grid.m(1);
37 m_y = obj.grid.m(2);
38 m_z = obj.grid.m(3);
39 obj.params = params; 33 obj.params = params;
40 34 dim = obj.grid.d;
41 switch operator 35 assert(dim == 3, 'Dimensions not correct')
42 case 'upwind'
43 ops_x = sbp.D1Upwind(m_x,xlim,order);
44 ops_y = sbp.D1Upwind(m_y,ylim,order);
45 ops_z = sbp.D1Upwind(m_z,zlim,order);
46 otherwise
47 ops_x = sbp.D2Standard(m_x,xlim,order);
48 ops_y = sbp.D2Standard(m_y,ylim,order);
49 ops_z = sbp.D2Standard(m_z,zlim,order);
50 end
51
52 x = obj.grid.x{1};
53 y = obj.grid.x{2};
54 z = obj.grid.x{3};
55 points = obj.grid.points(); 36 points = obj.grid.points();
56 X = points(:, 1); 37 m = obj.grid.m;
57 Y = points(:, 2); 38 for i = 1:dim
58 Z = points(:, 3); 39 ops{i} = opSet(m(i),obj.grid.lim{i},order);
59 40 x{i} = obj.grid.x{i};
60 obj.Yx = kr(y, ones(m_z,1)); 41 X{i} = points(:,i);
61 obj.Zx = kr(ones(m_y,1), z); 42 end
62 obj.Xy = kr(x, ones(m_z,1)); 43
63 obj.Zy = kr(ones(m_x,1), z); 44 obj.Yx = kr(x{2}, ones(m(3),1));
64 obj.Xz = kr(x, ones(m_y,1)); 45 obj.Zx = kr(ones(m(2),1), x{3});
65 obj.Yz = kr(ones(m_z,1), y); 46 obj.Xy = kr(x{1}, ones(m(3),1));
66 47 obj.Zy = kr(ones(m(1),1), x{3});
67 obj.Aevaluated = obj.evaluateCoefficientMatrix(A, X, Y, Z); 48 obj.Xz = kr(x{1}, ones(m(2),1));
68 obj.Bevaluated = obj.evaluateCoefficientMatrix(B, X, Y, Z); 49 obj.Yz = kr(ones(m(3),1), x{2});
69 obj.Cevaluated = obj.evaluateCoefficientMatrix(C, X, Y, Z); 50
70 obj.Eevaluated = obj.evaluateCoefficientMatrix(E, X, Y, Z); 51 obj.Aevaluated = obj.evaluateCoefficientMatrix(A, X{:});
52 obj.Bevaluated = obj.evaluateCoefficientMatrix(B, X{:});
53 obj.Cevaluated = obj.evaluateCoefficientMatrix(C, X{:});
54 obj.Eevaluated = obj.evaluateCoefficientMatrix(E, X{:});
71 55
72 obj.nEquations = length(A(obj.params,0,0,0)); 56 obj.nEquations = length(A(obj.params,0,0,0));
73 57
74 I_n = speye(obj.nEquations); 58 I_n = speye(obj.nEquations);
75 I_x = speye(m_x); 59 I_x = speye(m(1));
76 obj.I_x = I_x; 60 I_y = speye(m(2));
77 I_y = speye(m_y); 61 I_z = speye(m(3));
78 obj.I_y = I_y;
79 I_z = speye(m_z);
80 obj.I_z = I_z;
81 I_N = kr(I_n,I_x,I_y,I_z); 62 I_N = kr(I_n,I_x,I_y,I_z);
82 63
83 obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z); 64 obj.Hxi = kr(I_n, ops{1}.HI, I_y,I_z);
84 obj.Hx = ops_x.H; 65 obj.Hx = ops{1}.H;
85 obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z); 66 obj.Hyi = kr(I_n, I_x, ops{2}.HI,I_z);
86 obj.Hy = ops_y.H; 67 obj.Hy = ops{2}.H;
87 obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI); 68 obj.Hzi = kr(I_n, I_x,I_y, ops{3}.HI);
88 obj.Hz = ops_z.H; 69 obj.Hz = ops{3}.H;
89 70
90 obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z); 71 obj.e_w = kr(I_n, ops{1}.e_l, I_y,I_z);
91 obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z); 72 obj.e_e = kr(I_n, ops{1}.e_r, I_y,I_z);
92 obj.e_s = kr(I_n, I_x, ops_y.e_l,I_z); 73 obj.e_s = kr(I_n, I_x, ops{2}.e_l,I_z);
93 obj.e_n = kr(I_n, I_x, ops_y.e_r,I_z); 74 obj.e_n = kr(I_n, I_x, ops{2}.e_r,I_z);
94 obj.e_b = kr(I_n, I_x, I_y, ops_z.e_l); 75 obj.e_b = kr(I_n, I_x, I_y, ops{3}.e_l);
95 obj.e_t = kr(I_n, I_x, I_y, ops_z.e_r); 76 obj.e_t = kr(I_n, I_x, I_y, ops{3}.e_r);
96 77
97 obj.order = order; 78 obj.order = order;
98 79
99 switch operator 80 switch toString(opSet)
100 case 'upwind' 81 case 'sbp.D1Upwind'
101 alphaA = max(abs(eig(A(params, x(end), y(end), z(end))))); 82 alphaA = max(abs(eig(A(params, x{1}(end), x{2}(end), x{3}(end)))));
102 alphaB = max(abs(eig(B(params, x(end), y(end), z(end))))); 83 alphaB = max(abs(eig(B(params, x{1}(end), x{2}(end), x{3}(end)))));
103 alphaC = max(abs(eig(C(params, x(end), y(end), z(end))))); 84 alphaC = max(abs(eig(C(params, x{1}(end), x{2}(end), x{3}(end)))));
104 85
105 Ap = (obj.Aevaluated+alphaA*I_N)/2; 86 Ap = (obj.Aevaluated+alphaA*I_N)/2;
106 Am = (obj.Aevaluated-alphaA*I_N)/2; 87 Am = (obj.Aevaluated-alphaA*I_N)/2;
107 Dpx = kr(I_n, ops_x.Dp, I_y,I_z); 88 Dpx = kr(I_n, ops{1}.Dp, I_y,I_z);
108 Dmx = kr(I_n, ops_x.Dm, I_y,I_z); 89 Dmx = kr(I_n, ops{1}.Dm, I_y,I_z);
109 obj.D = -Am*Dpx; 90 obj.D = -Am*Dpx;
110 temp = Ap*Dmx; 91 temp = Ap*Dmx;
111 obj.D = obj.D-temp; 92 obj.D = obj.D-temp;
112 clear Ap Am Dpx Dmx 93 clear Ap Am Dpx Dmx
113 94
114 Bp = (obj.Bevaluated+alphaB*I_N)/2; 95 Bp = (obj.Bevaluated+alphaB*I_N)/2;
115 Bm = (obj.Bevaluated-alphaB*I_N)/2; 96 Bm = (obj.Bevaluated-alphaB*I_N)/2;
116 Dpy = kr(I_n, I_x, ops_y.Dp,I_z); 97 Dpy = kr(I_n, I_x, ops{2}.Dp,I_z);
117 Dmy = kr(I_n, I_x, ops_y.Dm,I_z); 98 Dmy = kr(I_n, I_x, ops{2}.Dm,I_z);
118 temp = Bm*Dpy; 99 temp = Bm*Dpy;
119 obj.D = obj.D-temp; 100 obj.D = obj.D-temp;
120 temp = Bp*Dmy; 101 temp = Bp*Dmy;
121 obj.D = obj.D-temp; 102 obj.D = obj.D-temp;
122 clear Bp Bm Dpy Dmy 103 clear Bp Bm Dpy Dmy
123 104
124 105
125 Cp = (obj.Cevaluated+alphaC*I_N)/2; 106 Cp = (obj.Cevaluated+alphaC*I_N)/2;
126 Cm = (obj.Cevaluated-alphaC*I_N)/2; 107 Cm = (obj.Cevaluated-alphaC*I_N)/2;
127 Dpz = kr(I_n, I_x, I_y,ops_z.Dp); 108 Dpz = kr(I_n, I_x, I_y,ops{3}.Dp);
128 Dmz = kr(I_n, I_x, I_y,ops_z.Dm); 109 Dmz = kr(I_n, I_x, I_y,ops{3}.Dm);
129 110
130 temp = Cm*Dpz; 111 temp = Cm*Dpz;
131 obj.D = obj.D-temp; 112 obj.D = obj.D-temp;
132 temp = Cp*Dmz; 113 temp = Cp*Dmz;
133 obj.D = obj.D-temp; 114 obj.D = obj.D-temp;
134 clear Cp Cm Dpz Dmz 115 clear Cp Cm Dpz Dmz
135 obj.D = obj.D-obj.Eevaluated; 116 obj.D = obj.D-obj.Eevaluated;
136 117
137 case 'standard' 118 case {'sbp.D2Standard', 'sbp.D2Variable', 'sbp.D4Standard', 'sbp.D4Variable'}
138 D1_x = kr(I_n, ops_x.D1, I_y,I_z); 119 D1_x = kr(I_n, ops{1}.D1, I_y,I_z);
139 D1_y = kr(I_n, I_x, ops_y.D1,I_z); 120 D1_y = kr(I_n, I_x, ops{2}.D1,I_z);
140 D1_z = kr(I_n, I_x, I_y,ops_z.D1); 121 D1_z = kr(I_n, I_x, I_y,ops{3}.D1);
141 obj.D = -obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated; 122 obj.D = -obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated;
142 otherwise 123 otherwise
143 error('Opperator not supported'); 124 error('Operator not supported');
144 end 125 end
145 end 126 end
146 127
147 % Closure functions return the operators applied to the own domain to close the boundary 128 % Closure functions return the operators applied to the own domain to close the boundary
148 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. 129 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.