Mercurial > repos > public > sbplib

comparison +scheme/Burgers1D.m @ 834:f1f0bf087e1c feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add support for artificial viscosity - Add option to use upwind operators for discretizing the scheme - Add option to use dissipation operators constructed from undivided differences when discretizing with narrow stencil operators
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 13 Sep 2018 14:21:47 +0200
parents 9f4c45a2d271
children 9e4e0576ca0f
comparison
equal deleted inserted replaced
833:9f4c45a2d271 834:f1f0bf087e1c
4 order % Order accuracy for the approximation 4 order % Order accuracy for the approximation
5 5
6 params 6 params
7 7
8 D % Non-stabalized scheme operator 8 D % Non-stabalized scheme operator
9 M % Derivative norm
10 H % Discrete norm 9 H % Discrete norm
11 Hi % Norm inverse 10 Hi % Norm inverse
12 e_l 11 e_l
13 e_r 12 e_r
14 d_l 13 d_l
15 d_r 14 d_r
16 end 15 end
17 16
18 methods 17 methods
19 function obj = Burgers1D(grid, pde_form, operator_type, order, params) 18 function obj = Burgers1D(grid, pde_form, operator_type, order, dissipation, params)
20 assert(grid.D == 1); 19 assert(grid.D == 1);
21 assert(grid.size() == length(params.eps)); 20 assert(grid.size() == length(params.eps));
22 m = grid.size(); 21 m = grid.size();
23 h = grid.scaling();
24 lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids. 22 lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids.
25 default_arg('pde_form','skew-symmetric'); 23 default_arg('pde_form','skew-symmetric');
26 default_arg('operator_type','narrow'); 24 default_arg('operator_type','narrow');
25 default_arg('dissipation','on');
27 26
28 switch operator_type 27 switch operator_type
29 case 'narrow' 28 case 'narrow'
30 ops = sbp.D4Variable(m, lim, order); 29 ops = sbp.D4Variable(m, lim, order);
31 D1 = ops.D1; 30 D1 = ops.D1;
32 D2 = ops.D2; 31 D2 = ops.D2;
32 DissipationOp = -1*sbp.dissipationOperator(m, order, ops.HI);
33 d_l = ops.d1_l';
34 d_r = ops.d1_r';
35 case 'upwind'
36 ops = sbp.D1Upwind(m, lim, order);
37 D1 = (ops.Dp + ops.Dm)/2;
38 %D2eps = @(eps) ops.Dp*diag(eps)*ops.Dm;
39 %D2eps = @(eps) ops.Dm*diag(eps)*ops.Dp;
40 D2 = @(eps) (ops.Dp*diag(eps)*ops.Dm + ops.Dm*diag(eps)*ops.Dp)/2;
41 DissipationOp = (ops.Dp-ops.Dm)/2;
42 d_l = D1;
43 d_r = D1;
33 otherwise 44 otherwise
34 error('Other operator types not yet supported', operator_type); 45 error('Other operator types not yet supported', operator_type);
35 end 46 end
36 47
48 %% TODO: Figure out how to evaluate viscosity as viscosity(v,t) here instead of parametrizing D on the viscosity.
37 switch pde_form 49 switch pde_form
38 case 'skew-symmetric' 50 case 'skew-symmetric'
39 D = @(v, viscosity) -1/3*v.*D1*v - 1/3*D1*v.^2 + D2(params.eps + viscosity)*v; 51 if (strcmp(dissipation,'on'))
52 D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v;
53 else
54 D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity))*v;
55 end
40 case 'conservative' 56 case 'conservative'
41 D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v; 57 if (strcmp(dissipation,'on'))
58 D = @(v, viscosity) -1/2*D1*v.^2 + (D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v;
59 else
60 D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v;
61 end
42 end 62 end
43 63
44 obj.grid = grid; 64 obj.grid = grid;
45 obj.order = order; 65 obj.order = order;
46 obj.params = params; 66 obj.params = params;
47 67
48 %% TODO: Figure out how to evaluate viscosity as viscosity(v,t) here instead of parametrizing D on the viscosity.
49 obj.D = D; 68 obj.D = D;
50 obj.M = ops.M;
51 obj.H = ops.H; 69 obj.H = ops.H;
52 obj.Hi = ops.HI; 70 obj.Hi = ops.HI;
53 obj.e_l = ops.e_l; 71 obj.e_l = ops.e_l;
54 obj.e_r = ops.e_r; 72 obj.e_r = ops.e_r;
55 obj.d_l = ops.d1_l; 73 obj.d_l = d_l;
56 obj.d_r = ops.d1_r; 74 obj.d_r = d_r;
57 end 75 end
58 76
59 % Closure functions return the opertors applied to the own doamin to close the boundary 77 % Closure functions return the opertors applied to the own doamin to close the boundary
60 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other domain. 78 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other domain.
61 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 79 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
69 [e, s, d] = obj.get_boundary_ops(boundary); 87 [e, s, d] = obj.get_boundary_ops(boundary);
70 switch type 88 switch type
71 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary 89 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary
72 case {'R','robin'} 90 case {'R','robin'}
73 p = s*obj.Hi*e; 91 p = s*obj.Hi*e;
74 closure = @(v, viscosity) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*(obj.params.eps + viscosity)*d'*v); 92 closure = @(v, viscosity) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*((obj.params.eps + viscosity).*d*v));
75 switch class(data) 93 switch class(data)
76 case 'double' 94 case 'double'
77 penalty = s*p*data; 95 penalty = s*p*data;
78 case 'function_handle' 96 case 'function_handle'
79 penalty = @(t) s*p*data(t); 97 penalty = @(t) s*p*data(t);