Mercurial > repos > public > sbplib

comparison +scheme/Burgers1D.m @ 853:cda996e64925 feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Minor renaming and clean up
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 12 Oct 2018 08:41:57 +0200
parents fbb8be3177c8
children
comparison
equal deleted inserted replaced
852:fbb8be3177c8 853:cda996e64925
18 function obj = Burgers1D(grid, pde_form, operator_type, order, dissipation, params) 18 function obj = Burgers1D(grid, pde_form, operator_type, order, dissipation, params)
19 assert(grid.D == 1); 19 assert(grid.D == 1);
20 assert(grid.size() == length(params.eps)); 20 assert(grid.size() == length(params.eps));
21 m = grid.size(); 21 m = grid.size();
22 lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids. 22 lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids.
23 default_arg('pde_form','skew-symmetric');
24 default_arg('operator_type','narrow');
25 default_arg('dissipation','on');
26
27 switch operator_type 23 switch operator_type
28 case 'narrow' 24 case 'narrow'
29 ops = sbp.D4Variable(m, lim, order); 25 ops = sbp.D4Variable(m, lim, order);
30 D1 = ops.D1; 26 D1 = ops.D1;
31 D2 = ops.D2; 27 D2 = ops.D2;
63 d_r = ops.e_r'*D1; 59 d_r = ops.e_r'*D1;
64 otherwise 60 otherwise
65 error('Other operator types not yet supported', operator_type); 61 error('Other operator types not yet supported', operator_type);
66 end 62 end
67 63
68 %% TODO: Figure out how to evaluate viscosity as viscosity(v,t) here instead of parametrizing D on the viscosity.
69 switch pde_form 64 switch pde_form
70 case 'skew-symmetric' 65 case 'skew-symmetric'
71 if (strcmp(dissipation,'on')) 66 if (strcmp(dissipation,'on'))
72 D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v; 67 D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v;
73 else 68 else
77 if (strcmp(dissipation,'on')) 72 if (strcmp(dissipation,'on'))
78 D = @(v, viscosity) -1/2*D1*v.^2 + (D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v; 73 D = @(v, viscosity) -1/2*D1*v.^2 + (D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v;
79 else 74 else
80 D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v; 75 D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v;
81 end 76 end
77 otherwise
78 error('Not supported', pde_form);
82 end 79 end
83 80
84 obj.grid = grid; 81 obj.grid = grid;
85 obj.order = order; 82 obj.order = order;
86 obj.params = params; 83 obj.params = params;
92 obj.e_r = ops.e_r; 89 obj.e_r = ops.e_r;
93 obj.d_l = d_l; 90 obj.d_l = d_l;
94 obj.d_r = d_r; 91 obj.d_r = d_r;
95 end 92 end
96 93
97 % Closure functions return the opertors applied to the own doamin to close the boundary 94 % Closure functions return the operators applied to the own doamin to close the boundary
98 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other domain. 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.
99 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 96 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
100 % type is a string specifying the type of boundary condition if there are several. 97 % type is a string specifying the type of boundary condition if there are several.
101 % data is a function returning the data that should be applied at the boundary. 98 % data is a function returning the data that should be applied at the boundary.
102 % neighbour_scheme is an instance of Scheme that should be interfaced to.
103 % neighbour_boundary is a string specifying which boundary to interface to.
104 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 99 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
105 default_arg('type','robin'); 100 default_arg('type','robin');
106 default_arg('data',0); 101 default_arg('data',0);
107 [e, s, d, i_b] = obj.get_boundary_ops(boundary); 102 [e, d, i_b, s] = obj.get_boundary_ops(boundary);
108 switch type 103 switch type
109 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary 104 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary
110 case {'R','robin'} 105 case {'R','robin'}
111 p = s*obj.Hi*e; 106 p = s*obj.Hi*e;
112 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)); 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));
116 case 'function_handle' 111 case 'function_handle'
117 penalty = @(t) s*p*data(t); 112 penalty = @(t) s*p*data(t);
118 otherwise 113 otherwise
119 error('Wierd data argument!') 114 error('Wierd data argument!')
120 end 115 end
121 % Unknown, boundary condition
122 otherwise 116 otherwise
123 error('No such boundary condition: type = %s',type); 117 error('No such boundary condition: type = %s',type);
124 end 118 end
125 end 119 end
126 120
127 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 121 % Ruturns the boundary ops, boundary index and sign for the boundary specified by the string boundary.
128 % The right boundary is considered the positive boundary 122 % The right boundary is considered the positive boundary
129 function [e, s, d, i_b] = get_boundary_ops(obj,boundary) 123 function [e, d, i_b, s] = get_boundary_ops(obj,boundary)
130 switch boundary 124 switch boundary
131 case 'l' 125 case 'l'
132 e = obj.e_l; 126 e = obj.e_l;
133 s = -1;
134 d = obj.d_l; 127 d = obj.d_l;
135 i_b = 1; 128 i_b = 1;
129 s = -1;
136 case 'r' 130 case 'r'
137 e = obj.e_r; 131 e = obj.e_r;
138 s = 1;
139 d = obj.d_r; 132 d = obj.d_r;
140 i_b = length(e); 133 i_b = length(e);
134 s = 1;
141 otherwise 135 otherwise
142 error('No such boundary: boundary = %s',boundary); 136 error('No such boundary: boundary = %s',boundary);
143 end 137 end
144 end 138 end
145 139