Mercurial > repos > public > sbplib

comparison +scheme/Burgers1d.m @ 1038:8537fdd6830a feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use opSets in Burgers1d. Note. Also changed from viscid formulation to inviscid formulation. Might be changed back later.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 18 Jan 2019 09:04:38 +0100
parents 2b9bdb22baec
children 0a5503a08a36
comparison
equal deleted inserted replaced
1037:2d7ba44340d0 1038:8537fdd6830a
1 classdef Burgers1d < scheme.Scheme 1 classdef Burgers1d < scheme.Scheme
2 properties 2 properties
3 grid % Physical grid 3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing
5 grid % Grid
4 order % Order accuracy for the approximation 6 order % Order accuracy for the approximation
7
8 H % Discrete norm
9 D
5 10
6 params 11 D1
7 12 Hi
8 D % Non-stabalized scheme operator
9 H % Discrete norm
10 Hi % Norm inverse
11 e_l 13 e_l
12 e_r 14 e_r
13 d_l
14 d_r
15 end 15 end
16 16
17 methods 17 methods
18 function obj = Burgers1d(grid, pde_form, operator_type, order, dissipation, params) 18 function obj = Burgers1d(g, order, pde_form, fluxSplitting, opSet)
19 assert(grid.D == 1); 19 default_arg('opSet',@sbp.D2Standard);
20 assert(grid.size() == length(params.eps)); 20 default_arg('fluxSplitting',@(v)max(abs(v)));
21 m = grid.size(); 21 assertType(g, 'grid.Cartesian');
22 lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids. 22
23 switch operator_type 23 m = g.size();
24 case 'narrow' 24 xl = g.getBoundary('l');
25 ops = sbp.D4Variable(m, lim, order); 25 xr = g.getBoundary('r');
26 D1 = ops.D1; 26 xlim = {xl, xr};
27 D2 = ops.D2; 27
28 if (strcmp(dissipation,'on')) 28 ops = opSet(m, xlim, order);
29 DissipationOp = -1*sbp.dissipationOperator(m, order, ops.HI); 29
30 end 30 if (isequal(opSet, @sbp.D1Upwind))
31 d_l = ops.d1_l'; 31 obj.D1 = (ops.Dp + ops.Dm)/2;
32 d_r = ops.d1_r'; 32 DissOp = (ops.Dm - ops.Dp)/2;
33 case 'upwind-' 33 else
34 ops = sbp.D1Upwind(m, lim, order); 34 obj.D1 = ops.D1;
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 35 end
63 36
64 switch pde_form 37 switch pde_form
65 case 'skew-symmetric' 38 case 'skew-symmetric'
66 if (strcmp(dissipation,'on')) 39 if (isequal(opSet, @sbp.D1Upwind))
67 D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v; 40 D = @(v) -(1/3*obj.D1*v.*v + (1/3*spdiag(v)*obj.D1 + fluxSplitting(v)*DissOp)*v);
68 else 41 else
69 D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity))*v; 42 D = @(v) -(1/3*obj.D1*v.*v + 1/3*spdiag(v)*obj.D1*v);
70 end 43 end
71 case 'conservative' 44 case 'conservative'
72 if (strcmp(dissipation,'on')) 45 if (isequal(opSet, @sbp.D1Upwind))
73 D = @(v, viscosity) -1/2*D1*v.^2 + (D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v; 46 D = @(v) -(1/2*obj.D1*v.*v + fluxSplitting(v)*DissOp*v);
74 else 47 else
75 D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v; 48 D = @(v) -(1/2*obj.D1*v.*v);
76 end 49 end
77 otherwise 50 otherwise
78 error('Not supported', pde_form); 51 error('Not supported', pde_form);
79 end 52 end
80 53
81 obj.grid = grid; 54 obj.grid = g;
82 obj.order = order;
83 obj.params = params;
84 55
85 obj.D = D; 56 obj.D = D;
86 obj.H = ops.H; 57 obj.H = ops.H;
87 obj.Hi = ops.HI; 58 obj.Hi = ops.HI;
59
88 obj.e_l = ops.e_l; 60 obj.e_l = ops.e_l;
89 obj.e_r = ops.e_r; 61 obj.e_r = ops.e_r;
90 obj.d_l = d_l; 62
91 obj.d_r = d_r; 63 obj.m = m;
64 obj.h = ops.h;
65 obj.order = order;
92 end 66 end
93 67
94 % Closure functions return the operators applied to the own doamin to close the boundary 68 % 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. 69 % 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'. 70 % 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. 71 % 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. 72 function [closure, penalty] = boundary_condition(obj, boundary, type)
99 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
100 default_arg('type','robin'); 73 default_arg('type','robin');
101 default_arg('data',0); 74 [e, index, s] = obj.get_boundary_ops(boundary);
102 [e, d, i_b, s] = obj.get_boundary_ops(boundary);
103 switch type 75 switch type
104 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary 76 % Stable dirhclet-like boundary conditions ((u+-abs(u))*u/3)) with +- at left/right boundary
105 case {'R','robin'} 77 case {'D', 'd', 'dirichlet', 'Dirichlet'}
106 p = s*obj.Hi*e; 78 tau = s*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)); 79 closure = @(v) obj.Hi*tau*(((v(index)-s*abs(v(index)))/3)*v(index));
108 switch class(data) 80 penalty = -obj.Hi*tau;
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 81 otherwise
117 error('No such boundary condition: type = %s',type); 82 error('No such boundary condition: type = %s',type);
118 end 83 end
119 end 84 end
120 85
121 % Ruturns the boundary ops, boundary index and sign for the boundary specified by the string boundary. 86 % Returns the boundary ops, boundary index and sign for the boundary specified by the string boundary.
122 % The right boundary is considered the positive boundary 87 % The right boundary is considered the positive boundary
123 function [e, d, i_b, s] = get_boundary_ops(obj,boundary) 88 function [e, index, s] = get_boundary_ops(obj,boundary)
124 switch boundary 89 switch boundary
125 case 'l' 90 case {'l','L','left','Left'}
126 e = obj.e_l; 91 e = obj.e_l;
127 d = obj.d_l; 92 index = 1;
128 i_b = 1;
129 s = -1; 93 s = -1;
130 case 'r' 94 case {'r','R','right','Right'}
131 e = obj.e_r; 95 e = obj.e_r;
132 d = obj.d_r; 96 index = length(e);
133 i_b = length(e);
134 s = 1; 97 s = 1;
135 otherwise 98 otherwise
136 error('No such boundary: boundary = %s',boundary); 99 error('No such boundary: boundary = %s',boundary);
137 end 100 end
138 end 101 end