Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst3dCurve.m @ 366:7ada2db63268 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Tried to speed up the matrix set-up a little
author Ylva Rydin <ylva.rydin@telia.com>
date Tue, 13 Dec 2016 09:28:33 +0100
parents 69b078cf8072
children 53abf04f5e4e
comparison
equal deleted inserted replaced
361:4a5d83324ec1 366:7ada2db63268
62 ops_eta = sbp.D1Upwind(m_eta,etalim,order); 62 ops_eta = sbp.D1Upwind(m_eta,etalim,order);
63 ops_zeta = sbp.D1Upwind(m_zeta,zetalim,order); 63 ops_zeta = sbp.D1Upwind(m_zeta,zetalim,order);
64 case 'standard' 64 case 'standard'
65 ops_xi = sbp.D2Standard(m_xi,xilim,order); 65 ops_xi = sbp.D2Standard(m_xi,xilim,order);
66 ops_eta = sbp.D2Standard(m_eta,etalim,order); 66 ops_eta = sbp.D2Standard(m_eta,etalim,order);
67 ops_zeta = sbp.D2Standard(m_zeta,zetalim,order); 67 ops_zeta = sbp.D2Standard(m_zeta,zetalim,order);
68 otherwise 68 otherwise
69 error('Operator not available') 69 error('Operator not available')
70 end 70 end
71 71
72 obj.xi = ops_xi.x; 72 obj.xi = ops_xi.x;
75 75
76 obj.Xi = kr(obj.xi,ones(m_eta,1),ones(m_zeta,1));%% Que pasa? 76 obj.Xi = kr(obj.xi,ones(m_eta,1),ones(m_zeta,1));%% Que pasa?
77 obj.Eta = kr(ones(m_xi,1),obj.eta,ones(m_zeta,1)); 77 obj.Eta = kr(ones(m_xi,1),obj.eta,ones(m_zeta,1));
78 obj.Zeta = kr(ones(m_xi,1),ones(m_eta,1),obj.zeta); 78 obj.Zeta = kr(ones(m_xi,1),ones(m_eta,1),obj.zeta);
79 79
80 obj.Eta_xi = kr(obj.eta,ones(m_xi,1));
81 obj.Zeta_xi = kr(ones(m_eta,1),obj.zeta);
82
83 obj.Xi_eta = kr(obj.xi,ones(m_zeta,1));
84 obj.Zeta_eta = kr(ones(m_xi,1),obj.zeta);
85
86 obj.Xi_zeta = kr(obj.xi,ones(m_eta,1));
87 obj.Eta_zeta = kr(ones(m_zeta,1),obj.eta);
88 80
89 [X,Y,Z] = ti.map(obj.Xi,obj.Eta,obj.Zeta); 81 [X,Y,Z] = ti.map(obj.Xi,obj.Eta,obj.Zeta);
90 obj.X = X; 82 obj.X = X;
91 obj.Y = Y; 83 obj.Y = Y;
92 obj.Z = Z; 84 obj.Z = Z;
95 I_xi = speye(m_xi); 87 I_xi = speye(m_xi);
96 obj.I_xi = I_xi; 88 obj.I_xi = I_xi;
97 I_eta = speye(m_eta); 89 I_eta = speye(m_eta);
98 obj.I_eta = I_eta; 90 obj.I_eta = I_eta;
99 I_zeta = speye(m_zeta); 91 I_zeta = speye(m_zeta);
100 obj.I_zeta = I_zeta; 92 obj.I_zeta = I_zeta;
101 93
102 I_N=kr(I_n,I_xi,I_eta,I_zeta); 94 I_N=kr(I_n,I_xi,I_eta,I_zeta);
103 95
104 O_xi = ones(m_xi,1); 96 O_xi = ones(m_xi,1);
105 O_eta = ones(m_eta,1); 97 O_eta = ones(m_eta,1);
106 O_zeta = ones(m_zeta,1); 98 O_zeta = ones(m_zeta,1);
107 99
108 obj.Hxii = kr(I_n, ops_xi.HI, I_eta,I_zeta); 100
109 obj.Hetai = kr(I_n, I_xi, ops_eta.HI,I_zeta);
110 obj.Hzetai = kr(I_n, I_xi,I_eta, ops_zeta.HI);
111 obj.Hxi = ops_xi.H; 101 obj.Hxi = ops_xi.H;
112 obj.Heta = ops_eta.H; 102 obj.Heta = ops_eta.H;
113 obj.Hzeta = ops_zeta.H; 103 obj.Hzeta = ops_zeta.H;
114 obj.h = [ops_xi.h ops_eta.h ops_zeta.h]; 104 obj.h = [ops_xi.h ops_eta.h ops_zeta.h];
115 105
116 switch operator 106 switch operator
117 case 'upwind' 107 case 'upwind'
118 D1_xi = kr((ops_xi.Dp+ops_xi.Dm)/2, I_eta,I_zeta); 108 D1_xi = kr((ops_xi.Dp+ops_xi.Dm)/2, I_eta,I_zeta);
119 D1_eta = kr(I_xi, (ops_eta.Dp+ops_eta.Dm)/2,I_zeta); 109 D1_eta = kr(I_xi, (ops_eta.Dp+ops_eta.Dm)/2,I_zeta);
120 D1_zeta = kr(I_xi, I_eta,(ops_zeta.Dp+ops_zeta.Dm)/2); 110 D1_zeta = kr(I_xi, I_eta,(ops_zeta.Dp+ops_zeta.Dm)/2);
121 otherwise 111 otherwise
122 D1_xi = kr(ops_xi.D1, I_eta,I_zeta); 112 D1_xi = kr(ops_xi.D1, I_eta,I_zeta);
123 D1_eta = kr(I_xi, ops_eta.D1,I_zeta); 113 D1_eta = kr(I_xi, ops_eta.D1,I_zeta);
124 D1_zeta = kr(I_xi, I_eta,ops_zeta.D1); 114 D1_zeta = kr(I_xi, I_eta,ops_zeta.D1);
125 end 115 end
126
127 obj.e_w = kr(I_n, ops_xi.e_l, I_eta,I_zeta);
128 obj.e_e = kr(I_n, ops_xi.e_r, I_eta,I_zeta);
129 obj.e_s = kr(I_n, I_xi, ops_eta.e_l,I_zeta);
130 obj.e_n = kr(I_n, I_xi, ops_eta.e_r,I_zeta);
131 obj.e_b = kr(I_n, I_xi, I_eta, ops_zeta.e_l);
132 obj.e_t = kr(I_n, I_xi, I_eta, ops_zeta.e_r);
133 116
134 obj.A = A; 117 obj.A = A;
135 obj.B = B; 118 obj.B = B;
136 obj.C = C; 119 obj.C = C;
137 120
142 obj.Y_eta = D1_eta*Y; 125 obj.Y_eta = D1_eta*Y;
143 obj.Y_zeta = D1_zeta*Y; 126 obj.Y_zeta = D1_zeta*Y;
144 obj.Z_xi = D1_xi*Z; 127 obj.Z_xi = D1_xi*Z;
145 obj.Z_eta = D1_eta*Z; 128 obj.Z_eta = D1_eta*Z;
146 obj.Z_zeta = D1_zeta*Z; 129 obj.Z_zeta = D1_zeta*Z;
147 130
148 D1_xi = kr(I_n,D1_xi); 131 obj.Ahat = @transform_coefficient_matrix;
149 D1_eta = kr(I_n,D1_eta); 132 obj.Bhat = @transform_coefficient_matrix;
150 D1_zeta = kr(I_n,D1_zeta); 133 obj.Chat = @transform_coefficient_matrix;
134 obj.E = @(obj,x,y,z,~,~,~,~,~,~)E(obj,x,y,z);
135
136 obj.Aevaluated = obj.evaluateCoefficientMatrix(obj.Ahat,obj.X, obj.Y,obj.Z, obj.X_eta,obj.X_zeta,obj.Y_eta,obj.Y_zeta,obj.Z_eta,obj.Z_zeta);
137 obj.Bevaluated = obj.evaluateCoefficientMatrix(obj.Bhat,obj.X, obj.Y,obj.Z, obj.X_zeta,obj.X_xi,obj.Y_zeta,obj.Y_xi,obj.Z_zeta,obj.Z_xi);
138 obj.Cevaluated = obj.evaluateCoefficientMatrix(obj.Chat,obj.X,obj.Y,obj.Z, obj.X_xi,obj.X_eta,obj.Y_xi,obj.Y_eta,obj.Z_xi,obj.Z_eta);
139
140 switch operator
141 case 'upwind'
142 clear D1_xi D1_eta D1_zeta
143 alphaA = max(eig(obj.Ahat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_eta(end),obj.X_zeta(end),obj.Y_eta(end),obj.Y_zeta(end),obj.Z_eta(end),obj.Z_zeta(end))));
144 alphaB = max(eig(obj.Bhat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_zeta(end),obj.X_xi(end),obj.Y_zeta(end),obj.Y_xi(end),obj.Z_zeta(end),obj.Z_xi(end))));
145 alphaC = max(eig(obj.Chat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_xi(end),obj.X_eta(end),obj.Y_xi(end),obj.Y_eta(end),obj.Z_xi(end),obj.Z_eta(end))));
146
147 Ap = (obj.Aevaluated+alphaA*I_N)/2;
148 Dmxi = kr(I_n, ops_xi.Dm, I_eta,I_zeta);
149 diffSum=-Ap*Dmxi;
150 clear Ap Dmxi
151
152 Am = (obj.Aevaluated-alphaA*I_N)/2;
153 obj.Aevaluated=[];
154 Dpxi = kr(I_n, ops_xi.Dp, I_eta,I_zeta);
155 temp=Am*Dpxi;
156 diffSum=diffSum-temp;
157 clear Am Dpxi
158
159 Bp = (obj.Bevaluated+alphaB*I_N)/2;
160 Dmeta = kr(I_n, I_xi, ops_eta.Dm,I_zeta);
161 temp=Bp*Dmeta;
162 diffSum=diffSum-temp;
163 clear Bp Dmeta
164
165 Bm = (obj.Bevaluated-alphaB*I_N)/2;
166 obj.Bevaluated=[];
167 Dpeta = kr(I_n, I_xi, ops_eta.Dp,I_zeta);
168 temp=Bm*Dpeta;
169 diffSum=diffSum-temp;
170 clear Bm Dpeta
171
172
173 Cp = (obj.Cevaluated+alphaC*I_N)/2;
174 Dmzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dm);
175 temp=Cp*Dmzeta;
176 diffSum=diffSum-temp;
177 clear Cp Dmzeta
178
179 Cm = (obj.Cevaluated-alphaC*I_N)/2;
180 clear I_N
181 obj.Cevaluated=[];
182 Dpzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dp);
183 temp=Cm*Dpzeta;
184 diffSum=diffSum-temp;
185 clear Cm Dpzeta temp
186
187
188 obj.J = obj.X_xi.*obj.Y_eta.*obj.Z_zeta...
189 +obj.X_zeta.*obj.Y_xi.*obj.Z_eta...
190 +obj.X_eta.*obj.Y_zeta.*obj.Z_xi...
191 -obj.X_xi.*obj.Y_zeta.*obj.Z_eta...
192 -obj.X_eta.*obj.Y_xi.*obj.Z_zeta...
193 -obj.X_zeta.*obj.Y_eta.*obj.Z_xi;
194
195 obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot));
196 obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]);
197
198 obj.D = obj.Ji*diffSum-obj.Eevaluated;
199 otherwise
200 D1_xi=kr(I_n,D1_xi);
201 D1_eta=kr(I_n,D1_eta);
202 D1_zeta=kr(I_n,D1_zeta);
203
204 obj.J = obj.X_xi.*obj.Y_eta.*obj.Z_zeta...
205 +obj.X_zeta.*obj.Y_xi.*obj.Z_eta...
206 +obj.X_eta.*obj.Y_zeta.*obj.Z_xi...
207 -obj.X_xi.*obj.Y_zeta.*obj.Z_eta...
208 -obj.X_eta.*obj.Y_xi.*obj.Z_zeta...
209 -obj.X_zeta.*obj.Y_eta.*obj.Z_xi;
210
211 obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot));
212 obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]);
213
214 obj.D = obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated;
215 end
216 obj.Hxii = kr(I_n, ops_xi.HI, I_eta,I_zeta);
217 obj.Hetai = kr(I_n, I_xi, ops_eta.HI,I_zeta);
218 obj.Hzetai = kr(I_n, I_xi,I_eta, ops_zeta.HI);
151 219
152 obj.index_w = (kr(ops_xi.e_l, O_eta,O_zeta)==1); 220 obj.index_w = (kr(ops_xi.e_l, O_eta,O_zeta)==1);
153 obj.index_e = (kr(ops_xi.e_r, O_eta,O_zeta)==1); 221 obj.index_e = (kr(ops_xi.e_r, O_eta,O_zeta)==1);
154 obj.index_s = (kr(O_xi, ops_eta.e_l,O_zeta)==1); 222 obj.index_s = (kr(O_xi, ops_eta.e_l,O_zeta)==1);
155 obj.index_n = (kr(O_xi, ops_eta.e_r,O_zeta)==1); 223 obj.index_n = (kr(O_xi, ops_eta.e_r,O_zeta)==1);
156 obj.index_b = (kr(O_xi, O_eta, ops_zeta.e_l)==1); 224 obj.index_b = (kr(O_xi, O_eta, ops_zeta.e_l)==1);
157 obj.index_t = (kr(O_xi, O_eta, ops_zeta.e_r)==1); 225 obj.index_t = (kr(O_xi, O_eta, ops_zeta.e_r)==1);
158 226
159 227 obj.e_w = kr(I_n, ops_xi.e_l, I_eta,I_zeta);
160 obj.Ahat = @transform_coefficient_matrix; 228 obj.e_e = kr(I_n, ops_xi.e_r, I_eta,I_zeta);
161 obj.Bhat = @transform_coefficient_matrix; 229 obj.e_s = kr(I_n, I_xi, ops_eta.e_l,I_zeta);
162 obj.Chat = @transform_coefficient_matrix; 230 obj.e_n = kr(I_n, I_xi, ops_eta.e_r,I_zeta);
163 obj.E = @(obj,x,y,z,~,~,~,~,~,~)E(obj,x,y,z); 231 obj.e_b = kr(I_n, I_xi, I_eta, ops_zeta.e_l);
164 232 obj.e_t = kr(I_n, I_xi, I_eta, ops_zeta.e_r);
165 obj.Aevaluated = obj.evaluateCoefficientMatrix(obj.Ahat,obj.X, obj.Y,obj.Z, obj.X_eta,obj.X_zeta,obj.Y_eta,obj.Y_zeta,obj.Z_eta,obj.Z_zeta); 233
166 obj.Bevaluated = obj.evaluateCoefficientMatrix(obj.Bhat,obj.X, obj.Y,obj.Z, obj.X_zeta,obj.X_xi,obj.Y_zeta,obj.Y_xi,obj.Z_zeta,obj.Z_xi); 234
167 obj.Cevaluated = obj.evaluateCoefficientMatrix(obj.Chat,obj.X,obj.Y,obj.Z, obj.X_xi,obj.X_eta,obj.Y_xi,obj.Y_eta,obj.Z_xi,obj.Z_eta); 235
168 obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); 236 obj.Eta_xi = kr(obj.eta,ones(m_xi,1));
169 237 obj.Zeta_xi = kr(ones(m_eta,1),obj.zeta);
170 obj.J = obj.X_xi.*obj.Y_eta.*obj.Z_zeta... 238 obj.Xi_eta = kr(obj.xi,ones(m_zeta,1));
171 +obj.X_zeta.*obj.Y_xi.*obj.Z_eta... 239 obj.Zeta_eta = kr(ones(m_xi,1),obj.zeta);
172 +obj.X_eta.*obj.Y_zeta.*obj.Z_xi... 240 obj.Xi_zeta = kr(obj.xi,ones(m_eta,1));
173 -obj.X_xi.*obj.Y_zeta.*obj.Z_eta... 241 obj.Eta_zeta = kr(ones(m_zeta,1),obj.eta);
174 -obj.X_eta.*obj.Y_xi.*obj.Z_zeta...
175 -obj.X_zeta.*obj.Y_eta.*obj.Z_xi;
176
177 obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot));
178
179
180 switch operator
181 case 'upwind'
182 alphaA = max(eig(obj.Ahat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_eta(end),obj.X_zeta(end),obj.Y_eta(end),obj.Y_zeta(end),obj.Z_eta(end),obj.Z_zeta(end))));
183 alphaB = max(eig(obj.Bhat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_zeta(end),obj.X_xi(end),obj.Y_zeta(end),obj.Y_xi(end),obj.Z_zeta(end),obj.Z_xi(end))));
184 alphaC = max(eig(obj.Chat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_xi(end),obj.X_eta(end),obj.Y_xi(end),obj.Y_eta(end),obj.Z_xi(end),obj.Z_eta(end))));
185
186 Ap = (obj.Aevaluated+alphaA*I_N)/2;
187 Dmxi = kr(I_n, ops_xi.Dm, I_eta,I_zeta);
188 diffSum=-Ap*Dmxi;
189 clear Ap Dmxi
190
191 Am = (obj.Aevaluated-alphaA*I_N)/2;
192 clear obj.Aevaluated
193 Dpxi = kr(I_n, ops_xi.Dp, I_eta,I_zeta);
194 temp=Am*Dpxi;
195 diffSum=diffSum-temp;
196 clear Am Dpxi
197
198 Bp = (obj.Bevaluated+alphaB*I_N)/2;
199 Dmeta = kr(I_n, I_xi, ops_eta.Dm,I_zeta);
200 temp=Bp*Dmeta;
201 diffSum=diffSum-temp;
202 clear Bp Dmeta
203
204 Bm = (obj.Bevaluated-alphaB*I_N)/2;
205 clear obj.Bevaluated
206 Dpeta = kr(I_n, I_xi, ops_eta.Dp,I_zeta);
207 temp=Bm*Dpeta;
208 diffSum=diffSum-temp;
209
210
211 Cp = (obj.Cevaluated+alphaC*I_N)/2;
212 Dmzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dm);
213 temp=Cp*Dmzeta;
214 diffSum=diffSum-temp;
215 clear Cp Dmzeta
216
217 Cm = (obj.Cevaluated-alphaC*I_N)/2;
218 clear obj.Cevaluated
219 Dpzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dp);
220 temp=Cm*Dpzeta;
221 diffSum=diffSum-temp;
222
223 obj.D = obj.Ji*diffSum-obj.Eevaluated;
224 otherwise
225 obj.D = obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated;
226 end
227
228 242
229 end 243 end
230 244
231 function [ret] = transform_coefficient_matrix(obj,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2) 245 function [ret] = transform_coefficient_matrix(obj,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2)
232 ret = obj.A(obj,x,y,z).*(y_1.*z_2-z_1.*y_2); 246 ret = obj.A(obj,x,y,z).*(y_1.*z_2-z_1.*y_2);
279 matVec = mat; 293 matVec = mat;
280 [rows,cols] = size(matVec); 294 [rows,cols] = size(matVec);
281 side = max(length(X),length(Y)); 295 side = max(length(X),length(Y));
282 cols = cols/side; 296 cols = cols/side;
283 end 297 end
298 matVec(abs(matVec)<10^(-10))=0;
284 ret = cell(rows,cols); 299 ret = cell(rows,cols);
285 300
286 301
287 for ii=1:rows 302 for ii=1:rows
288 for jj=1:cols 303 for jj=1:cols
289 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); 304 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side));
290 end 305 end
291 end 306 end
292 307
293 ret = cell2mat(ret); 308 ret = cell2mat(ret);
309
294 end 310 end
295 311
296 312
297 function [BM] = boundary_matrices(obj,boundary) 313 function [BM] = boundary_matrices(obj,boundary)
298 params = obj.params; 314 params = obj.params;