Mercurial > repos > public > sbplib

comparison +scheme/Burgers2d.m @ 1149:1fe48cbd379a feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Change Burgers2d to inviscid formulation. Rewrite to use opSets and fix the implementation of the Dirichlet conditions.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 24 Jan 2019 09:05:44 +0100
parents c322a9454d75
children 3108963cc42c 0c906f7ab8bf
comparison
equal deleted inserted replaced
1148:0a5503a08a36 1149:1fe48cbd379a
1 classdef Burgers2D < scheme.Scheme 1 classdef Burgers2d < scheme.Scheme
2 properties 2 properties
3 grid % Physical grid 3 grid % Physical grid
4 order % Order accuracy for the approximation 4 order % Order accuracy for the approximation
5 5
6 D % Non-stabilized scheme operator 6 D % Non-stabilized scheme operator
7 H % Discrete norm 7 H % Discrete norm
8 H_inv % Norm inverse 8 H_x, H_y % Norms in the x and y directions
9 halfnorm_inv % Cell array halfnorm operators 9 Hi % Kroneckered norm inverse
10 e_l % Cell array of left boundary operators 10 % Boundary operators
11 e_r % Cell array of right boundary operators 11 e_w, e_e, e_s, e_n
12 d_l % Cell array of left boundary derivative operators
13 d_r % Cell array of right boundary derivative operators
14 end 12 end
15 13
16 methods 14 methods
17 function obj = Burgers2D(g, operator_type, order, dissipation) 15 function obj = Burgers2d(g, order, pde_form, fluxSplitting, opSet)
18 if ~isa(g, 'grid.Cartesian') || g.D() ~= 2 16 default_arg('opSet',@sbp.D2Standard);
19 error('Grid must be 2d cartesian'); 17 default_arg('fluxSplitting',{@(v)max(abs(v)),@(v)max(abs(v))});
18 assertType(g, 'grid.Cartesian');
19
20 m = g.size();
21 m_x = m(1);
22 m_y = m(2);
23 m_tot = g.N();
24
25 xlim = {g.x{1}(1), g.x{1}(end)};
26 ylim = {g.x{2}(1), g.x{2}(end)};
27 obj.grid = g;
28
29 % Operator sets
30 ops_x = opSet(m_x, xlim, order);
31 ops_y = opSet(m_y, ylim, order);
32 Ix = speye(m_x);
33 Iy = speye(m_y);
34
35 % Norms
36 Hx = ops_x.H;
37 Hy = ops_y.H;
38 Hxi = ops_x.HI;
39 Hyi = ops_y.HI;
40
41 obj.H_x = Hx;
42 obj.H_y = Hy;
43 obj.H = kron(Hx,Hy);
44 obj.Hi = kron(Hxi,Hyi);
45
46 % Derivatives
47 if (isequal(opSet,@sbp.D1Upwind))
48 Dx = kron((ops_x.Dp + ops_x.Dm)/2,Iy);
49 Dy = kron(Ix,(ops_y.Dp + ops_y.Dm)/2);
50 DissOpx = kron((ops_x.Dm - ops_x.Dp)/2,Iy);
51 DissOpy = kron(Ix,(ops_y.Dm - ops_y.Dp)/2);
52 D1 = Dx + Dy;
53 switch pde_form
54 case 'skew-symmetric'
55 switch length(fluxSplitting)
56 case 1
57 DissOp = DissOpx + DissOpy;
58 obj.D = @(v) -(1/3*D1*v.*v + (1/3*spdiag(v)*D1 + fluxSplitting{1}(v)*DissOp)*v);
59 case 2
60 obj.D = @(v) -(1/3*D1*v.*v + (1/3*spdiag(v)*D1 + fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v);
61 end
62 case 'conservative'
63 switch length(fluxSplitting)
64 case 1
65 DissOp = DissOpx + DissOpy;
66 obj.D = @(v) -(1/2*D1*v.*v + fluxSplitting{1}(v)*DissOp*v);
67 case 2
68 obj.D = @(v) -(1/2*D1*v.*v + (fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v);
69 end
70
71 end
72 else
73 Dx = kron(ops_x.D1,Iy);
74 Dy = kron(Ix,ops_y.D1);
75 D1 = Dx + Dy;
76 switch pde_form
77 case 'skew-symmetric'
78 obj.D = @(v) -(1/3*D1*v.*v + 1/3*spdiag(v)*D1*v);
79 case 'conservative'
80 obj.D = @(v) -1/2*D1*v.*v;
81 end
20 end 82 end
21 83
22 obj.grid = g; 84 % Boundary operators
85 obj.e_w = kr(ops_x.e_l, Iy);
86 obj.e_e = kr(ops_x.e_r, Iy);
87 obj.e_s = kr(Ix, ops_y.e_l);
88 obj.e_n = kr(Ix, ops_y.e_r);
89
23 obj.order = order; 90 obj.order = order;
24
25 % Create cell array of 1D operators. For example D1_1d{1} = D1_x, D1_1d{2} = D1_y.
26 [Dp_1d, Dm_1d, H_1d, H_inv_1d, d_l_1d, d_r_1d, e_l_1d, e_r_1d, I, DissipationOp_1d] = ...
27 obj.assemble1DOperators(g, operator_type, order, dissipation);
28
29 %% 2D-operators
30 % D1
31 D1_1d{1} = (Dp_1d{1} + Dp_1d{1})/2;
32 D1_1d{2} = (Dp_1d{2} + Dp_1d{2})/2;
33 D1_2d = obj.extendOperatorTo2D(D1_1d, I);
34 D1 = D1_2d{1} + D1_2d{2};
35 % D2
36
37 Dp_2d = obj.extendOperatorTo2D(Dp_1d, I);
38 Dm_2d = obj.extendOperatorTo2D(Dm_1d, I);
39 D2 = @(viscosity) Dm_2d{1}*spdiag(viscosity)*Dp_2d{1} + Dm_2d{2}*spdiag(viscosity)*Dp_2d{2};
40 % m = g.size();
41 % ind = grid.funcToMatrix(g, 1:g.N());
42 % for i = 1:g.D()
43 % D2_2d{i} = sparse(zeros(g.N()));
44 % end
45 % % x-direction
46 % for i = 1:m(2)
47 % p = ind(:,i);
48 % D2_2d{1}(p,p) = @(viscosity) D2_1d{1}(viscosity(p));
49 % end
50 % % y-direction
51 % for i = 1:m(1)
52 % p = ind(i,:);
53 % D2_2d{2}(p,p) = @(viscosity) D2_1d{2}(viscosity(p));
54 % end
55 % D2 = D2_2d{1} + D2_2d{2};
56
57 obj.d_l = obj.extendOperatorTo2D(d_l_1d, I);
58 obj.d_r = obj.extendOperatorTo2D(d_r_1d, I);
59 obj.e_l = obj.extendOperatorTo2D(e_l_1d, I);
60 obj.e_r = obj.extendOperatorTo2D(e_r_1d, I);
61 obj.H = kron(H_1d{1},H_1d{2});
62 obj.H_inv = kron(H_inv_1d{1},H_inv_1d{2});
63 obj.halfnorm_inv = obj.extendOperatorTo2D(H_inv_1d, I);
64
65 % Dissipation operator
66 switch dissipation
67 case 'on'
68 DissOp_2d = obj.extendOperatorTo2D(DissipationOp_1d, I);
69 DissOp = DissOp_2d{1} + DissOp_2d{2};
70 obj.D = @(v, viscosity) -1/2*D1*v.^2 + (D2(viscosity) + max(abs(v))*DissOp)*v;
71 case 'off'
72 obj.D = @(v, viscosity) -1/2*D1*v.^2 + D2(viscosity)*v;
73 end
74 end 91 end
75 92
76 % Closure functions return the operators applied to the own doamin to close the boundary 93 % Closure functions return the operators applied to the own doamin to close the boundary
77 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain. 94 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain.
78 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 95 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
79 % type is a string specifying the type of boundary condition if there are several. 96 % type is a string specifying the type of boundary condition if there are several.
80 % data is a function returning the data that should be applied at the boundary. 97 function [closure, penalty] = boundary_condition(obj,boundary,type)
81 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 98 default_arg('type','dirichlet');
82 default_arg('type','robin'); 99 [e, H_b, index, s] = obj.get_boundary_ops(boundary);
83 default_arg('data',0);
84 [e, d, halfnorm_inv, i_b, s] = obj.get_boundary_ops(boundary);
85 switch type 100 switch type
86 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary 101 % Stable dirchlet-like boundary conditions (u+-abs(u))*u/3
87 case {'R','robin'} 102 % with +- at left/right boundaries in each coordinate direction
88 p = s*halfnorm_inv*e; 103 case {'D', 'd', 'dirichlet', 'Dirichlet'}
89 closure = @(v, viscosity) p*(((v(i_b)-s*abs(v(i_b)))/3).*(v(i_b)) - ((viscosity(i_b)).*d*v)); 104
90 switch class(data) 105 magnitude = 2/3;
91 case 'double' 106 tau = @(v) s*magnitude*obj.Hi*e*H_b*spdiag((v(index)-s*abs(v(index)))/2);
92 penalty = s*p*data; 107 closure = @(v) tau(v)*v(index);
93 case 'function_handle' 108 penalty = @(v) -tau(v);
94 penalty = @(t) s*p*data(t); 109
95 otherwise 110
96 error('Wierd data argument!') 111 % tau = s*e*H_b;
97 end 112 % closure = @(v) (obj.Hi*tau)*(((v(index)-s*abs(v(index)))/3).*v(index));
113 % penalty = -obj.Hi*tau;
98 otherwise 114 otherwise
99 error('No such boundary condition: type = %s',type); 115 error('No such boundary condition: type = %s',type);
100 end 116 end
117
118
101 end 119 end
102 120
103 % Ruturns the boundary ops, half-norm, boundary indices and sign for the boundary specified by the string boundary. 121 % Ruturns the boundary ops, half-norm, boundary indices and sign for the boundary specified by the string boundary.
104 % The right boundary for each coordinate direction is considered the positive boundary 122 % The right boundary for each coordinate direction is considered the positive boundary
105 function [e, d, halfnorm_inv, ind_boundary, s] = get_boundary_ops(obj,boundary) 123 function [e, H_b, index, s] = get_boundary_ops(obj, boundary)
106 ind = grid.funcToMatrix(obj.grid, 1:obj.grid.N()); 124 ind = grid.funcToMatrix(obj.grid, 1:obj.grid.N());
107 switch boundary 125 switch boundary
108 case 'w' 126 case {'w', 'W', 'west', 'West'}
109 e = obj.e_l{1}; 127 e = obj.e_w;
110 d = obj.d_l{1}; 128 H_b = obj.H_y;
111 halfnorm_inv = obj.halfnorm_inv{1}; 129 index = ind(1,:);
112 ind_boundary = ind(1,:);
113 s = -1; 130 s = -1;
114 case 'e' 131 case {'e', 'E', 'east', 'East'}
115 e = obj.e_r{1}; 132 e = obj.e_e;
116 d = obj.d_r{1}; 133 H_b = obj.H_y;
117 halfnorm_inv = obj.halfnorm_inv{1}; 134 index = ind(end,:);
118
119 ind_boundary = ind(end,:);
120 s = 1; 135 s = 1;
121 case 's' 136 case {'s', 'S', 'south', 'South'}
122 e = obj.e_l{2}; 137 e = obj.e_s;
123 d = obj.d_l{2}; 138 H_b = obj.H_x;
124 halfnorm_inv = obj.halfnorm_inv{2}; 139 index = ind(:,1);
125 ind_boundary = ind(:,1);
126 s = -1; 140 s = -1;
127 case 'n' 141 case {'n', 'N', 'north', 'North'}
128 e = obj.e_r{2}; 142 e = obj.e_n;
129 d = obj.d_r{2}; 143 H_b = obj.H_x;
130 halfnorm_inv = obj.halfnorm_inv{2}; 144 index = ind(:,end);
131 ind_boundary = ind(:,end);
132 s = 1; 145 s = 1;
133 otherwise 146 otherwise
134 error('No such boundary: boundary = %s',boundary); 147 error('No such boundary: boundary = %s',boundary);
135 end 148 end
136 end 149 end
141 154
142 function N = size(obj) 155 function N = size(obj)
143 N = obj.grid.m; 156 N = obj.grid.m;
144 end 157 end
145 end 158 end
146
147 methods(Static)
148 function [Dp, Dm, H, Hi, d_l, d_r, e_l, e_r, I, DissipationOp] = assemble1DOperators(g, operator_type, order, dissipation)
149 dim = g.D();
150 I = cell(dim,1);
151 D1 = cell(dim,1);
152 D2 = cell(dim,1);
153 H = cell(dim,1);
154 Hi = cell(dim,1);
155 e_l = cell(dim,1);
156 e_r = cell(dim,1);
157 d1_l = cell(dim,1);
158 d1_r = cell(dim,1);
159 DissipationOp = cell(dim,1);
160 for i = 1:dim
161 switch operator_type
162 % case 'narrow'
163 % ops = sbp.D4Variable(g.m(i), g.lim{i}, order);
164 % D1{i} = ops.D1;
165 % D2{i} = ops.D2;
166 % d_l{i} = ops.d1_l';
167 % d_r{i} = ops.d1_r';
168 % if (strcmp(dissipation,'on'))
169 % DissipationOp{i} = -1*sbp.dissipationOperator(g.m(i), order, ops.HI);
170 % end
171 % case 'upwind-'
172 % ops = sbp.D1Upwind(g.m(i), g.lim{i}, order);
173 % D1{i} = (ops.Dp + ops.Dm)/2;
174 % D2{i} = @(viscosity) ops.Dp*spdiag(viscosity)*ops.Dm;
175 % d_l{i} = ops.e_l'*ops.Dm;
176 % d_r{i} = ops.e_r'*ops.Dm;
177 % if (strcmp(dissipation,'on'))
178 % DissipationOp{i} = (ops.Dp-ops.Dm)/2;
179 % end
180 case 'upwind+'
181 ops = sbp.D1Upwind(g.m(i), g.lim{i}, order);
182 Dp{i} = ops.Dp;
183 Dm{i} = ops.Dm;
184 % D1{i} = (ops.Dp + ops.Dm)/2;
185 % D2{i} = @(viscosity) ops.Dm*spdiag(viscosity)*ops.Dp;
186 d_l{i} = ops.e_l'*ops.Dp;
187 d_r{i} = ops.e_r'*ops.Dp;
188 if (strcmp(dissipation,'on'))
189 DissipationOp{i} = (ops.Dp-ops.Dm)/2;
190 end
191 % case 'upwind+-'
192 % ops = sbp.D1Upwind(g.m(i), g.lim{i}, order);
193 % D1{i} = (ops.Dp + ops.Dm)/2;
194 % D2{i} = @(viscosity) (ops.Dp*spdiag(viscosity)*ops.Dm + ops.Dm*spdiag(viscosity)*ops.Dp)/2;
195 % d_l{i} = ops.e_l'*D1;
196 % d_r{i} = ops.e_r'*D1;
197 % if (strcmp(dissipation,'on'))
198 % DissipationOp{i} = (ops.Dp-ops.Dm)/2;
199 % end
200 otherwise
201 error('Other operator types not yet supported', operator_type);
202 end
203 H{i} = ops.H;
204 Hi{i} = ops.HI;
205 e_l{i} = ops.e_l;
206 e_r{i} = ops.e_r;
207 I{i} = speye(g.m(i));
208 end
209 end
210 function op_2d = extendOperatorTo2D(op, I)
211 op_2d{1} = kr(op{1}, I{2});
212 op_2d{2} = kr(I{1}, op{2});
213 end
214 end
215 end 159 end