Mercurial > repos > public > sbplib

comparison +scheme/Burgers2D.m @ 1008:a6f34de60044 feature/burgers2d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at implementing Burgers in 2D with RV-stabilization
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 12 Oct 2018 08:54:39 +0200
parents
children
comparison
equal deleted inserted replaced
854:18162a0a5bb5 1008:a6f34de60044
1 classdef Burgers2D < scheme.Scheme
2 properties
3 grid % Physical grid
4 order % Order accuracy for the approximation
5
6 D % Non-stabilized scheme operator
7 H % Discrete norm
8 H_inv % Norm inverse
9 halfnorm_inv % Cell array halfnorm operators
10 e_l % Cell array of left boundary operators
11 e_r % Cell array of right boundary operators
12 d_l % Cell array of left boundary derivative operators
13 d_r % Cell array of right boundary derivative operators
14 end
15
16 methods
17 function obj = Burgers2D(g, operator_type, order, dissipation)
18 if ~isa(g, 'grid.Cartesian') || g.D() ~= 2
19 error('Grid must be 2d cartesian');
20 end
21
22 obj.grid = g;
23 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
75
76 % 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.
78 % 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.
80 % data is a function returning the data that should be applied at the boundary.
81 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
82 default_arg('type','robin');
83 default_arg('data',0);
84 [e, d, halfnorm_inv, i_b, s] = obj.get_boundary_ops(boundary);
85 switch type
86 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary
87 case {'R','robin'}
88 p = s*halfnorm_inv*e;
89 closure = @(v, viscosity) p*(((v(i_b)-s*abs(v(i_b)))/3).*(v(i_b)) - ((viscosity(i_b)).*d*v));
90 switch class(data)
91 case 'double'
92 penalty = s*p*data;
93 case 'function_handle'
94 penalty = @(t) s*p*data(t);
95 otherwise
96 error('Wierd data argument!')
97 end
98 otherwise
99 error('No such boundary condition: type = %s',type);
100 end
101 end
102
103 % 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
105 function [e, d, halfnorm_inv, ind_boundary, s] = get_boundary_ops(obj,boundary)
106 ind = grid.funcToMatrix(obj.grid, 1:obj.grid.N());
107 switch boundary
108 case 'w'
109 e = obj.e_l{1};
110 d = obj.d_l{1};
111 halfnorm_inv = obj.halfnorm_inv{1};
112 ind_boundary = ind(1,:);
113 s = -1;
114 case 'e'
115 e = obj.e_r{1};
116 d = obj.d_r{1};
117 halfnorm_inv = obj.halfnorm_inv{1};
118
119 ind_boundary = ind(end,:);
120 s = 1;
121 case 's'
122 e = obj.e_l{2};
123 d = obj.d_l{2};
124 halfnorm_inv = obj.halfnorm_inv{2};
125 ind_boundary = ind(:,1);
126 s = -1;
127 case 'n'
128 e = obj.e_r{2};
129 d = obj.d_r{2};
130 halfnorm_inv = obj.halfnorm_inv{2};
131 ind_boundary = ind(:,end);
132 s = 1;
133 otherwise
134 error('No such boundary: boundary = %s',boundary);
135 end
136 end
137
138 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
139 error('An interface function does not exist yet');
140 end
141
142 function N = size(obj)
143 N = obj.grid.m;
144 end
145 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