Mercurial > repos > public > sbplib

comparison +scheme/Burgers1d.m @ 1035:2b9bdb22baec feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Forgot to add .m as file extension
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 17 Jan 2019 11:06:13 +0100
parents +scheme/Burgers1d@2676ad79f994
children 8537fdd6830a
comparison
equal deleted inserted replaced
1034:2676ad79f994 1035:2b9bdb22baec
1 classdef Burgers1d < scheme.Scheme
2 properties
3 grid % Physical grid
4 order % Order accuracy for the approximation
5
6 params
7
8 D % Non-stabalized scheme operator
9 H % Discrete norm
10 Hi % Norm inverse
11 e_l
12 e_r
13 d_l
14 d_r
15 end
16
17 methods
18 function obj = Burgers1d(grid, pde_form, operator_type, order, dissipation, params)
19 assert(grid.D == 1);
20 assert(grid.size() == length(params.eps));
21 m = grid.size();
22 lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids.
23 switch operator_type
24 case 'narrow'
25 ops = sbp.D4Variable(m, lim, order);
26 D1 = ops.D1;
27 D2 = ops.D2;
28 if (strcmp(dissipation,'on'))
29 DissipationOp = -1*sbp.dissipationOperator(m, order, ops.HI);
30 end
31 d_l = ops.d1_l';
32 d_r = ops.d1_r';
33 case 'upwind-'
34 ops = sbp.D1Upwind(m, lim, order);
35 D1 = (ops.Dp + ops.Dm)/2;
36 D2 = @(eps) ops.Dp*spdiag(eps)*ops.Dm;
37 if (strcmp(dissipation,'on'))
38 DissipationOp = (ops.Dp-ops.Dm)/2;
39 end
40 d_l = ops.e_l'*ops.Dm;
41 d_r = ops.e_r'*ops.Dm;
42 case 'upwind+'
43 ops = sbp.D1Upwind(m, lim, order);
44 D1 = (ops.Dp + ops.Dm)/2;
45 D2 = @(eps) ops.Dm*spdiag(eps)*ops.Dp;
46 if (strcmp(dissipation,'on'))
47 DissipationOp = (ops.Dp-ops.Dm)/2;
48 end
49 d_l = ops.e_l'*ops.Dp;
50 d_r = ops.e_r'*ops.Dp;
51 case 'upwind+-'
52 ops = sbp.D1Upwind(m, lim, order);
53 D1 = (ops.Dp + ops.Dm)/2;
54 D2 = @(eps) (ops.Dp*spdiag(eps)*ops.Dm + ops.Dm*spdiag(eps)*ops.Dp)/2;
55 if (strcmp(dissipation,'on'))
56 DissipationOp = (ops.Dp-ops.Dm)/2;
57 end
58 d_l = ops.e_l'*D1;
59 d_r = ops.e_r'*D1;
60 otherwise
61 error('Other operator types not yet supported', operator_type);
62 end
63
64 switch pde_form
65 case 'skew-symmetric'
66 if (strcmp(dissipation,'on'))
67 D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v;
68 else
69 D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity))*v;
70 end
71 case 'conservative'
72 if (strcmp(dissipation,'on'))
73 D = @(v, viscosity) -1/2*D1*v.^2 + (D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v;
74 else
75 D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v;
76 end
77 otherwise
78 error('Not supported', pde_form);
79 end
80
81 obj.grid = grid;
82 obj.order = order;
83 obj.params = params;
84
85 obj.D = D;
86 obj.H = ops.H;
87 obj.Hi = ops.HI;
88 obj.e_l = ops.e_l;
89 obj.e_r = ops.e_r;
90 obj.d_l = d_l;
91 obj.d_r = d_r;
92 end
93
94 % Closure functions return the operators applied to the own doamin to close the boundary
95 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain.
96 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
97 % type is a string specifying the type of boundary condition if there are several.
98 % data is a function returning the data that should be applied at the boundary.
99 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
100 default_arg('type','robin');
101 default_arg('data',0);
102 [e, d, i_b, s] = obj.get_boundary_ops(boundary);
103 switch type
104 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary
105 case {'R','robin'}
106 p = s*obj.Hi*e;
107 closure = @(v, viscosity) p*(((v(i_b)-s*abs(v(i_b)))/3)*(v(i_b)) - ((obj.params.eps(i_b) + viscosity(i_b))*d*v));
108 switch class(data)
109 case 'double'
110 penalty = s*p*data;
111 case 'function_handle'
112 penalty = @(t) s*p*data(t);
113 otherwise
114 error('Wierd data argument!')
115 end
116 otherwise
117 error('No such boundary condition: type = %s',type);
118 end
119 end
120
121 % Ruturns the boundary ops, boundary index and sign for the boundary specified by the string boundary.
122 % The right boundary is considered the positive boundary
123 function [e, d, i_b, s] = get_boundary_ops(obj,boundary)
124 switch boundary
125 case 'l'
126 e = obj.e_l;
127 d = obj.d_l;
128 i_b = 1;
129 s = -1;
130 case 'r'
131 e = obj.e_r;
132 d = obj.d_r;
133 i_b = length(e);
134 s = 1;
135 otherwise
136 error('No such boundary: boundary = %s',boundary);
137 end
138 end
139
140 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
141 error('An interface function does not exist yet');
142 end
143
144 function N = size(obj)
145 N = obj.grid.m;
146 end
147 end
148 end