Mercurial > repos > public > sbplib

comparison +scheme/Utux2D.m @ 595:2a2f34778ded feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make Utux2D work
author Martin Almquist <malmquist@stanford.edu>
date Tue, 26 Sep 2017 13:26:25 -0700
parents 39554f2de783
children 0f9d20dbb7ce
comparison
equal deleted inserted replaced
594:a2ddaccf5fd1 595:2a2f34778ded
7 v0 % Initial data 7 v0 % Initial data
8 8
9 a % Wave speed a = [a1, a2]; 9 a % Wave speed a = [a1, a2];
10 10
11 H % Discrete norm 11 H % Discrete norm
12 Hi, Hx, Hy, Hxi, Hyi 12 H_x, H_y % Norms in the x and y directions
13 Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms
13 14
14 % Derivatives 15 % Derivatives
15 Dx, Dy 16 Dx, Dy
16 17
17 % Boundary operators 18 % Boundary operators
32 m = g.size(); 33 m = g.size();
33 m_x = m(1); 34 m_x = m(1);
34 m_y = m(2); 35 m_y = m(2);
35 m_tot = g.N(); 36 m_tot = g.N();
36 37
37 xlim = g.x{1}; 38 xlim = {g.x{1}(1), g.x{1}(end)};
38 ylim = g.x{2}; 39 ylim = {g.x{2}(1), g.x{2}(end)};
39 obj.grid = g; 40 obj.grid = g;
40 41
41 % Operator sets 42 % Operator sets
42 ops_x = opSet(m_x, xlim, order); 43 ops_x = opSet(m_x, xlim, order);
43 ops_y = opSet(m_y, ylim, order); 44 ops_y = opSet(m_y, ylim, order);
47 % Norms 48 % Norms
48 Hx = ops_x.H; 49 Hx = ops_x.H;
49 Hy = ops_y.H; 50 Hy = ops_y.H;
50 Hxi = ops_x.HI; 51 Hxi = ops_x.HI;
51 Hyi = ops_y.HI; 52 Hyi = ops_y.HI;
53
54 obj.H_x = Hx;
55 obj.H_y = Hy;
52 obj.H = kron(Hx,Hy); 56 obj.H = kron(Hx,Hy);
53 obj.Hi = kron(Hxi,Hyi); 57 obj.Hi = kron(Hxi,Hyi);
54 obj.Hx = kron(Hx,Iy); 58 obj.Hx = kron(Hx,Iy);
55 obj.Hy = kron(Ix,Hy); 59 obj.Hy = kron(Ix,Hy);
56 obj.Hxi = kron(Hxi,Iy); 60 obj.Hxi = kron(Hxi,Iy);
69 obj.e_n = kr(Ix, ops_y.e_r); 73 obj.e_n = kr(Ix, ops_y.e_r);
70 74
71 obj.m = m; 75 obj.m = m;
72 obj.h = [ops_x.h ops_y.h]; 76 obj.h = [ops_x.h ops_y.h];
73 obj.order = order; 77 obj.order = order;
74 78 obj.a = a;
75 obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy); 79 obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy);
76 80
77 end 81 end
78 % Closure functions return the opertors applied to the own domain to close the boundary 82 % Closure functions return the opertors applied to the own domain to close the boundary
79 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 83 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
86 default_arg('type','dirichlet'); 90 default_arg('type','dirichlet');
87 91
88 sigma = -1; % Scalar penalty parameter 92 sigma = -1; % Scalar penalty parameter
89 switch boundary 93 switch boundary
90 case {'w','W','west','West'} 94 case {'w','W','west','West'}
91 tau = sigma*obj.a(1)*obj.e_w*obj.Hy; 95 tau = sigma*obj.a(1)*obj.e_w*obj.H_y;
92 closure = obj.Hi*tau*obj.e_w'; 96 closure = obj.Hi*tau*obj.e_w';
93 97
94 case {'s','S','south','South'} 98 case {'s','S','south','South'}
95 tau = sigma*pbj.a(2)*obj.e_s*obj.Hx; 99 tau = sigma*obj.a(2)*obj.e_s*obj.H_x;
96 closure = obj.Hi*tau*obj.e_s'; 100 closure = obj.Hi*tau*obj.e_s';
97 end 101 end
98 penalty = -obj.Hi*tau; 102 penalty = -obj.Hi*tau;
99 103
100 end 104 end
117 sigma_ds = -1; %"Downstream" penalty 121 sigma_ds = -1; %"Downstream" penalty
118 sigma_us = 0; %"Upstream" penalty 122 sigma_us = 0; %"Upstream" penalty
119 123
120 switch boundary 124 switch boundary
121 case {'w','W','west','West'} 125 case {'w','W','west','West'}
122 tau = sigma_ds*obj.a(1)*obj.e_w*obj.Hy; 126 tau = sigma_ds*obj.a(1)*obj.e_w*obj.H_y;
123 closure = obj.Hi*tau*obj.e_w'; 127 closure = obj.Hi*tau*obj.e_w';
124 case {'e','E','east','East'} 128 case {'e','E','east','East'}
125 tau = sigma_us*obj.a(1)*obj.e_e*obj.Hy; 129 tau = sigma_us*obj.a(1)*obj.e_e*obj.H_y;
126 closure = obj.Hi*tau*obj.e_e'; 130 closure = obj.Hi*tau*obj.e_e';
127 case {'s','S','south','South'} 131 case {'s','S','south','South'}
128 tau = sigma_ds*obj.a(2)*obj.e_s*obj.Hx; 132 tau = sigma_ds*obj.a(2)*obj.e_s*obj.H_x;
129 closure = obj.Hi*tau*obj.e_s'; 133 closure = obj.Hi*tau*obj.e_s';
130 case {'n','N','north','North'} 134 case {'n','N','north','North'}
131 tau = sigma_us*obj.a(2)*obj.e_n*obj.Hx; 135 tau = sigma_us*obj.a(2)*obj.e_n*obj.H_x;
132 closure = obj.Hi*tau*obj.e_n'; 136 closure = obj.Hi*tau*obj.e_n';
133 end 137 end
134 penalty = -obj.Hi*tau*e_neighbour'; 138 penalty = -obj.Hi*tau*e_neighbour';
135 139
136 end 140 end