Mercurial > repos > public > sbplib

comparison +scheme/Burgers1D.m @ 848:c8ea9bbdc62c feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Increase efficiency of computing closures
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 20 Sep 2018 19:02:01 +0200
parents 1e057b0f2fed
children ab2e5a24ddde
comparison
equal deleted inserted replaced
847:1c6f1595bb94 848:c8ea9bbdc62c
86 % neighbour_scheme is an instance of Scheme that should be interfaced to. 86 % neighbour_scheme is an instance of Scheme that should be interfaced to.
87 % neighbour_boundary is a string specifying which boundary to interface to. 87 % neighbour_boundary is a string specifying which boundary to interface to.
88 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 88 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
89 default_arg('type','robin'); 89 default_arg('type','robin');
90 default_arg('data',0); 90 default_arg('data',0);
91 [e, s, d] = obj.get_boundary_ops(boundary); 91 [e, s, d, i_b] = obj.get_boundary_ops(boundary);
92 switch type 92 switch type
93 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary 93 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary
94 case {'R','robin'} 94 case {'R','robin'}
95 p = s*obj.Hi*e; 95 p = s*obj.Hi*e;
96 closure = @(v, viscosity) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*((obj.params.eps + viscosity).*d*v)); 96 closure = @(v, viscosity) p*(((v(i_b)-s*abs(v(i_b)))/3)*(v(i_b)) - e'*((obj.params.eps(i_b) + viscosity(i_b))*d*v));
97 switch class(data) 97 switch class(data)
98 case 'double' 98 case 'double'
99 penalty = s*p*data; 99 penalty = s*p*data;
100 case 'function_handle' 100 case 'function_handle'
101 penalty = @(t) s*p*data(t); 101 penalty = @(t) s*p*data(t);
108 end 108 end
109 end 109 end
110 110
111 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 111 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
112 % The right boundary is considered the positive boundary 112 % The right boundary is considered the positive boundary
113 function [e, s, d] = get_boundary_ops(obj,boundary) 113 function [e, s, d, i_b] = get_boundary_ops(obj,boundary)
114 switch boundary 114 switch boundary
115 case 'l' 115 case 'l'
116 e = obj.e_l; 116 e = obj.e_l;
117 s = -1; 117 s = -1;
118 d = obj.d_l; 118 d = obj.d_l;
119 i_b = 1;
119 case 'r' 120 case 'r'
120 e = obj.e_r; 121 e = obj.e_r;
121 s = 1; 122 s = 1;
122 d = obj.d_r; 123 d = obj.d_r;
124 i_b = length(e);
123 otherwise 125 otherwise
124 error('No such boundary: boundary = %s',boundary); 126 error('No such boundary: boundary = %s',boundary);
125 end 127 end
126 end 128 end
127 129