Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2dCurve.m @ 491:26125cfefe11 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Schrodinger 2d runs without exploding with fixed coordinates
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 10 Feb 2017 14:29:53 +0100
parents +grid/Schrodinger2dCurve.m@b13d44271ead
children 6b95a894cbd7
comparison
equal deleted inserted replaced
490:b13d44271ead 491:26125cfefe11
1 classdef Schrodinger2dCurve < scheme.Scheme
2 properties
3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing
5
6 grid
7 xm, ym
8
9 order % Order accuracy for the approximation
10
11 D % non-stabalized scheme operator
12 M % Derivative norm
13 H % Discrete norm
14 Hi
15 H_u, H_v % Norms in the x and y directions
16 Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
17 Hi_u, Hi_v
18 Hiu, Hiv
19 D1_v, D1_u
20 D2_v, D2_u
21 Du, Dv
22
23
24 e_w, e_e, e_s, e_n
25 du_w, dv_w
26 du_e, dv_e
27 du_s, dv_s
28 du_n, dv_n
29 g_1, g_2
30 c
31 a11, a12, a22
32 m_tot, m_u, m_v
33 p,p_tau
34 Ji
35 end
36
37 methods
38 function obj = Schrodinger2dCurve(g ,order, opSet,p,p_tau)
39 default_arg('opSet',@sbp.D2Variable);
40 default_arg('c', 1);
41
42 obj.p=p;
43 obj.p_tau=p_tau;
44 obj.c=1;
45
46 m = g.size();
47 obj.m_u = m(1);
48 obj.m_v = m(2);
49 obj.m_tot = g.N();
50
51 h = g.scaling();
52 h_u = h(1);
53 h_v = h(2);
54
55 % Operators
56 ops_u = opSet(obj.m_u, {0, 1}, order);
57 ops_v = opSet(obj.m_v, {0, 1}, order);
58
59 I_u = speye(obj.m_u);
60 I_v = speye(obj.m_v);
61
62 obj.D1_u = ops_u.D1;
63 obj.D2_u = ops_u.D2;
64
65 H_u = ops_u.H;
66 Hi_u = ops_u.HI;
67 e_l_u = ops_u.e_l;
68 e_r_u = ops_u.e_r;
69 d1_l_u = ops_u.d1_l;
70 d1_r_u = ops_u.d1_r;
71
72 obj.D1_v = ops_v.D1;
73 obj.D2_v = ops_v.D2;
74 H_v = ops_v.H;
75 Hi_v = ops_v.HI;
76 e_l_v = ops_v.e_l;
77 e_r_v = ops_v.e_r;
78 d1_l_v = ops_v.d1_l;
79 d1_r_v = ops_v.d1_r;
80
81 obj.Du = kr(obj.D1_u,I_v);
82 obj.Dv = kr(I_u,obj.D1_v);
83
84 obj.H = kr(H_u,H_v);
85 obj.Hi = kr(Hi_u,Hi_v);
86 obj.Hu = kr(H_u,I_v);
87 obj.Hv = kr(I_u,H_v);
88 obj.Hiu = kr(Hi_u,I_v);
89 obj.Hiv = kr(I_u,Hi_v);
90
91 obj.e_w = kr(e_l_u,I_v);
92 obj.e_e = kr(e_r_u,I_v);
93 obj.e_s = kr(I_u,e_l_v);
94 obj.e_n = kr(I_u,e_r_v);
95 obj.du_w = kr(d1_l_u,I_v);
96 obj.dv_w = (obj.e_w'*obj.Dv)';
97 obj.du_e = kr(d1_r_u,I_v);
98 obj.dv_e = (obj.e_e'*obj.Dv)';
99 obj.du_s = (obj.e_s'*obj.Du)';
100 obj.dv_s = kr(I_u,d1_l_v);
101 obj.du_n = (obj.e_n'*obj.Du)';
102 obj.dv_n = kr(I_u,d1_r_v);
103
104 % obj.x_u = x_u;
105 % obj.x_v = x_v;
106 % obj.y_u = y_u;
107 % obj.y_v = y_v;
108
109 obj.m = m;
110 obj.h = [h_u h_v];
111 obj.order = order;
112 obj.grid = g;
113
114
115 end
116
117
118 function [D ]= d_fun(obj,t)
119 % Metric derivatives
120 ti = parametrization.Ti.points(obj.p.p1(t),obj.p.p2(t),obj.p.p3(t),obj.p.p4(t));
121 ti_tau = parametrization.Ti.points(obj.p_tau.p1(t),obj.p_tau.p2(t),obj.p_tau.p3(t),obj.p_tau.p4(t));
122
123 lcoords=points(obj.grid);
124 [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2));
125 [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2));
126 x = reshape(obj.xm,obj.m_tot,1);
127 y = reshape(obj.ym,obj.m_tot,1);
128
129 x_tau = reshape(x_tau,obj.m_tot,1);
130 y_tau = reshape(y_tau,obj.m_tot,1);
131
132 x_u = obj.Du*x;
133 x_v = obj.Dv*x;
134 y_u = obj.Du*y;
135 y_v = obj.Dv*y;
136
137 J = x_u.*y_v - x_v.*y_u;
138 a11 = 1./J.* (x_v.^2 + y_v.^2);
139 a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
140 a22 = 1./J .* (x_u.^2 + y_u.^2);
141
142 obj.a11=a11;
143 obj.a12=a12;
144 obj.a22=a22;
145
146 % Assemble full operators
147 L_12 = spdiags(a12, 0, obj.m_tot, obj.m_tot);
148 Duv = obj.Du*L_12*obj.Dv;
149 Dvu = obj.Dv*L_12*obj.Du;
150
151 Duu = sparse(obj.m_tot);
152 Dvv = sparse(obj.m_tot);
153 ind = grid.funcToMatrix(obj.grid, 1:obj.m_tot);
154
155
156 for i = 1:obj.m_v
157 D = obj.D2_u(a11(ind(:,i)));
158 p = ind(:,i);
159 Duu(p,p) = D;
160 end
161
162 for i = 1:obj.m_u
163 D = obj.D2_v(a22(ind(i,:)));
164 p = ind(i,:);
165 Dvv(p,p) = D;
166 end
167
168 Ji = spdiags(1./J, 0, obj.m_tot, obj.m_tot);
169 obj.Ji=Ji;
170 obj.g_1 = x_tau.*y_v-y_tau.*x_v;
171 obj.g_2 = -x_tau.*y_u + y_tau.*x_u;
172
173 %Add the flux splitting
174 D = Ji*(obj.g_1.*obj.Du + obj.g_2.*obj.Dv + 1i*obj.c^2*(Duu + Duv + Dvu + Dvv)); %% g_1' och g_2'?
175
176 % obj.gamm_u = h_u*ops_u.borrowing.M.d1;
177 % obj.gamm_v = h_v*ops_v.borrowing.M.d1;
178
179 end
180
181 % Closure functions return the opertors applied to the own doamin to close the boundary
182 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
183 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
184 % type is a string specifying the type of boundary condition if there are several.
185 % data is a function returning the data that should be applied at the boundary.
186 % neighbour_scheme is an instance of Scheme that should be interfaced to.
187 % neighbour_boundary is a string specifying which boundary to interface to.
188 function [closure, penalty] = boundary_condition(obj, boundary)
189 [e, d_n, d_t, coeff_t, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary);
190
191 a_t = spdiag(coeff_t);
192 F = (s * d_n' + s * a_t*d_t')';
193 tau1 = 1;
194 a = spdiag(g);
195 tau2 = (-1*s*a - abs(a))/4;
196
197 penalty_parameter_1 = 1*1i*halfnorm_inv_n*halfnorm_inv_t*e*F'*halfnorm_t*e;
198 penalty_parameter_2 = halfnorm_inv_n*e*tau2;
199
200 closure = obj.Ji*obj.c^2 * penalty_parameter_1*e' +obj.Ji* penalty_parameter_2*e';
201 % penalty = -obj.Ji*obj.c^2 * penalty_parameter_2;
202 penalty=0;
203
204 end
205
206 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
207 end
208
209 function [e, d_n, d_t, coeff_t, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g, I] = get_boundary_ops(obj, boundary)
210
211 % gridMatrix = zeros(obj.m(2),obj.m(1));
212 % gridMatrix(:) = 1:numel(gridMatrix);
213
214 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
215
216 switch boundary
217 case 'w'
218 e = obj.e_w;
219 d_n = obj.du_w;
220 d_t = obj.dv_w;
221 s = -1;
222
223 I = ind(1,:);
224 coeff_t = obj.a12(I);
225 g=obj.g_1(I);
226 case 'e'
227 e = obj.e_e;
228 d_n = obj.du_e;
229 d_t = obj.dv_e;
230 s = 1;
231
232 I = ind(end,:);
233 coeff_t = obj.a12(I);
234 g=obj.g_1(I);
235 case 's'
236 e = obj.e_s;
237 d_n = obj.dv_s;
238 d_t = obj.du_s;
239 s = -1;
240
241 I = ind(:,1)';
242 coeff_t = obj.a12(I);
243 g=obj.g_2(I);
244 case 'n'
245 e = obj.e_n;
246 d_n = obj.dv_n;
247 d_t = obj.du_n;
248 s = 1;
249
250 I = ind(:,end)';
251 coeff_t = obj.a12(I);
252 g=obj.g_2(I);
253 otherwise
254 error('No such boundary: boundary = %s',boundary);
255 end
256
257 switch boundary
258 case {'w','e'}
259 halfnorm_inv_n = obj.Hiu;
260 halfnorm_inv_t = obj.Hiv;
261 halfnorm_t = obj.Hv;
262
263 case {'s','n'}
264 halfnorm_inv_n = obj.Hiv;
265 halfnorm_inv_t = obj.Hiu;
266 halfnorm_t = obj.Hu;
267 end
268 end
269
270 function N = size(obj)
271 N = prod(obj.m);
272 end
273
274
275 end
276 end