Mercurial > repos > public > sbplib

comparison +scheme/Utux2d.m @ 1025:ac80bedc8df7 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up of Utux2d
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 07 Jan 2019 16:26:05 +0100
parents 234c1c02ea39
children 8a9393084b30
comparison
equal deleted inserted replaced
1024:81d02e2d5c23 1025:ac80bedc8df7
2 properties 2 properties
3 m % Number of points in each direction, possibly a vector 3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing 4 h % Grid spacing
5 grid % Grid 5 grid % Grid
6 order % Order accuracy for the approximation 6 order % Order accuracy for the approximation
7 v0 % Initial data
8 7
9 a % Wave speed a = [a1, a2]; 8 a % Wave speed a = [a1, a2];
10 % Can either be a constant vector or a cell array of function handles. 9 % Can either be a constant vector or a cell array of function handles.
11 10
12 H % Discrete norm 11 H % Discrete norm
13 H_x, H_y % Norms in the x and y directions 12 H_x, H_y % Norms in the x and y directions
14 Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms 13 Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms
15 14
16 % Derivatives 15 % Derivatives
17 Dx, Dy 16 Dx, Dy
18
19 % Dissipation Operators
20 DissOpx, DissOpy
21 17
22 % Boundary operators 18 % Boundary operators
23 e_w, e_e, e_s, e_n 19 e_w, e_e, e_s, e_n
24 20
25 D % Total discrete operator 21 D % Total discrete operator
74 70
75 % Derivatives 71 % Derivatives
76 if (isequal(opSet,@sbp.D1Upwind)) 72 if (isequal(opSet,@sbp.D1Upwind))
77 Dx = (ops_x.Dp + ops_x.Dm)/2; 73 Dx = (ops_x.Dp + ops_x.Dm)/2;
78 Dy = (ops_y.Dp + ops_y.Dm)/2; 74 Dy = (ops_y.Dp + ops_y.Dm)/2;
75 obj.Dx = kron(Dx,Iy);
76 obj.Dy = kron(Ix,Dy);
79 DissOpx = (ops_x.Dm - ops_x.Dp)/2; 77 DissOpx = (ops_x.Dm - ops_x.Dp)/2;
80 DissOpy = (ops_y.Dm - ops_y.Dp)/2; 78 DissOpy = (ops_y.Dm - ops_y.Dp)/2;
81 obj.Dx = kron(Dx,Iy); 79 DissOpx = kron(DissOpx,Iy);
82 obj.Dy = kron(Ix,Dy); 80 DissOpy = kron(Ix,DissOpy);
83 obj.DissOpx = kron(DissOpx,Iy); 81
84 obj.DissOpy = kron(Ix,DissOpy); 82 obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy + fluxSplitting{1}*DissOpx + fluxSplitting{2}*DissOpy);
85 else 83 else
86 Dx = ops_x.D1; 84 Dx = ops_x.D1;
87 Dy = ops_y.D1; 85 Dy = ops_y.D1;
88 obj.Dx = kron(Dx,Iy); 86 obj.Dx = kron(Dx,Iy);
89 obj.Dy = kron(Ix,Dy); 87 obj.Dy = kron(Ix,Dy);
88
89 obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy);
90 end 90 end
91 91
92 % Boundary operators 92 % Boundary operators
93 obj.e_w = kr(ops_x.e_l, Iy); 93 obj.e_w = kr(ops_x.e_l, Iy);
94 obj.e_e = kr(ops_x.e_r, Iy); 94 obj.e_e = kr(ops_x.e_r, Iy);
97 97
98 obj.m = m; 98 obj.m = m;
99 obj.h = [ops_x.h ops_y.h]; 99 obj.h = [ops_x.h ops_y.h];
100 obj.order = order; 100 obj.order = order;
101 obj.a = a; 101 obj.a = a;
102
103 if (isequal(opSet,@sbp.D1Upwind))
104 obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy + fluxSplitting{1}*obj.DissOpx + fluxSplitting{2}*obj.DissOpy);
105 else
106 obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy);
107 end
108 end 102 end
109 % Closure functions return the opertors applied to the own domain to close the boundary 103 % Closure functions return the opertors applied to the own domain to close the boundary
110 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 104 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
111 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 105 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
112 % type is a string specifying the type of boundary condition if there are several. %TBD Remove type here? Only dirichlet applicable? 106 % type is a string specifying the type of boundary condition if there are several. %TBD Remove type here? Only dirichlet applicable?
113 % data is a function returning the data that should be applied at the boundary. 107 % data is a function returning the data that should be applied at the boundary.
114 % neighbour_scheme is an instance of Scheme that should be interfaced to. 108 % neighbour_scheme is an instance of Scheme that should be interfaced to.
115 % neighbour_boundary is a string specifying which boundary to interface to. 109 % neighbour_boundary is a string specifying which boundary to interface to.
116 function [closure, penalty] = boundary_condition(obj,boundary,type) 110 function [closure, penalty] = boundary_condition(obj,boundary,type)
117 default_arg('type','dirichlet'); 111 default_arg('type','dirichlet');
118 sigma_left = -1; % Scalar penalty parameter for left boundaries 112 sigma_left = -1; % Scalar penalty parameter for left boundaries (West/South)
119 sigma_right = 1; % Scalar penalty parameter for right boundaries 113 sigma_right = 1; % Scalar penalty parameter for right boundaries (East/North)
120 switch boundary 114 switch boundary
115 % Can only specify boundary condition where there is inflow
116 % Extract the postivie resp. negative part of a, for the left
117 % resp. right boundaries, and set other values of a to zero.
118 % Then the closure will effectively only contribute to inflow boundaries
121 case {'w','W','west','West'} 119 case {'w','W','west','West'}
122 % Can only specify boundary condition where there is inflow 120 a_inflow = obj.a{1};
123 % Extract positive part of a{1} 121 a_inflow(a_inflow < 0) = 0;
124 a_pos = obj.a{1}; 122 tau = sigma_left*a_inflow*obj.e_w*obj.H_y;
125 a_pos(a_pos < 0) = 0;
126 tau = sigma_left*a_pos*obj.e_w*obj.H_y;
127 closure = obj.Hi*tau*obj.e_w'; 123 closure = obj.Hi*tau*obj.e_w';
124 case {'e','E','east','East'}
125 a_inflow = obj.a{1};
126 a_inflow(a_inflow > 0) = 0;
127 tau = sigma_right*a_inflow*obj.e_e*obj.H_y;
128 closure = obj.Hi*tau*obj.e_e';
128 case {'s','S','south','South'} 129 case {'s','S','south','South'}
129 % Can only specify boundary condition where there is inflow 130 a_inflow = obj.a{2};
130 % Extract positive part of a{2} 131 a_inflow(a_inflow < 0) = 0;
131 a_pos = obj.a{2}; 132 tau = sigma_left*a_inflow*obj.e_s*obj.H_x;
132 a_pos(a_pos < 0) = 0;
133 tau = sigma_left*a_pos*obj.e_s*obj.H_x;
134 closure = obj.Hi*tau*obj.e_s'; 133 closure = obj.Hi*tau*obj.e_s';
135 case {'e','E','east','East'}
136 % Can only specify boundary condition where there is inflow
137 % Extract negative part of a{1}
138 a_neg = obj.a{1};
139 a_neg(a_neg > 0) = 0;
140 tau = sigma_right*a_neg*obj.e_e*obj.H_y;
141 closure = obj.Hi*tau*obj.e_e';
142 case {'n','N','north','North'} 134 case {'n','N','north','North'}
143 % Can only specify boundary condition where there is inflow 135 a_inflow = obj.a{2};
144 % Extract negative part of a{2} 136 a_inflow(a_inflow > 0) = 0;
145 a_neg = obj.a{2}; 137 tau = sigma_right*a_inflow*obj.e_n*obj.H_x;
146 a_neg(a_neg > 0) = 0;
147 tau = sigma_right*a_neg*obj.e_n*obj.H_x;
148 closure = obj.Hi*tau*obj.e_n'; 138 closure = obj.Hi*tau*obj.e_n';
149 end 139 end
150 penalty = -obj.Hi*tau; 140 penalty = -obj.Hi*tau;
151 end 141 end
152 142