Mercurial > repos > public > sbplib

comparison +scheme/Burgers1D.m @ 851:ab2e5a24ddde feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
- Fix bug when constructing closure for narrow stencils - Update the residual outside of the RK time steps. At least for the inner convergence rate, updating the residual inside does not seem to be required.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 21 Sep 2018 15:33:15 +0200
parents c8ea9bbdc62c
children fbb8be3177c8
comparison
equal deleted inserted replaced
850:4e5e53d6336c 851:ab2e5a24ddde
41 %D2 = @(eps) ops.Dm*diag(eps)*ops.Dp; 41 %D2 = @(eps) ops.Dm*diag(eps)*ops.Dp;
42 %D2 = @(eps) (ops.Dp*diag(eps)*ops.Dm + ops.Dm*diag(eps)*ops.Dp)/2; 42 %D2 = @(eps) (ops.Dp*diag(eps)*ops.Dm + ops.Dm*diag(eps)*ops.Dp)/2;
43 if (strcmp(dissipation,'on')) 43 if (strcmp(dissipation,'on'))
44 DissipationOp = (ops.Dp-ops.Dm)/2; 44 DissipationOp = (ops.Dp-ops.Dm)/2;
45 end 45 end
46 d_l = D1; 46 d_l = ops.e_l'*D1;
47 d_r = D1; 47 d_r = ops.e_r'*D1;
48 otherwise 48 otherwise
49 error('Other operator types not yet supported', operator_type); 49 error('Other operator types not yet supported', operator_type);
50 end 50 end
51 51
52 %% TODO: Figure out how to evaluate viscosity as viscosity(v,t) here instead of parametrizing D on the viscosity. 52 %% TODO: Figure out how to evaluate viscosity as viscosity(v,t) here instead of parametrizing D on the viscosity.
91 [e, s, d, i_b] = 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*(((v(i_b)-s*abs(v(i_b)))/3)*(v(i_b)) - e'*((obj.params.eps(i_b) + viscosity(i_b))*d*v)); 96 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));
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);