Mercurial > repos > public > sbplib
changeset 834:f1f0bf087e1c feature/burgers1d
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 | 008496ca38f3 |
| files | +sbp/dissipationOperator.m +scheme/Burgers1D.m |
| diffstat | 2 files changed, 97 insertions(+), 11 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/dissipationOperator.m Thu Sep 13 14:21:47 2018 +0200 @@ -0,0 +1,68 @@ +%% Function that constructs artificial dissipation operators using undivided differences +function D = dissipationOperator(m, order, Hinv, scaling) + % TBD: Add or remove D_2 and Dp/Dm? + % d2=[1 2 1]; + % D_2=(diag(ones(m-1,1),-1)-2*diag(ones(m,1),0)+ ... + % diag(ones(m-1,1),1)); + % D_2(1,1:3)=[d2];D_2(m,m-2:m)=[d2]; + % %Dm + % DD_m=(diag(ones(m-1,1),+1)-diag(ones(m,1),0)); + % DD_m(m,m-1:m)=[-1 1]; + % %Dp + % DD_p=(-diag(ones(m-1,1),-1)+diag(ones(m,1),0)); + % DD_p(1,1:2)=[-1 1]; + + switch order + case 1 + DD_1=(diag(ones(m-1,1),+1)-diag(ones(m,1),0)); + DD_1(m,m-1:m)=[-1 1]; + D = DD_2'*DD_2; + case 2 + dd2=0*[1 -2 1]; + DD_2=(diag(ones(m-1,1),-1)-2*diag(ones(m,1),0)+ ... + diag(ones(m-1,1),1)); + DD_2(1,1:3)=[dd2];DD_2(m,m-2:m)=[dd2]; + D = DD_2'*DD_2; + case 3 + d3=0*[-1 3 -3 1]; + DD_3=(-diag(ones(m-2,1),-2)+3*diag(ones(m-1,1),-1)-3*diag(ones(m,1),0)+ ... + diag(ones(m-1,1),1)); + DD_3(1:2,1:4)=[d3;d3]; + DD_3(m,m-3:m)=[d3]; + D = DD_3'*DD_3; + case 4 + default_arg('scaling', 1/12); + d4=0*[1 -4 6 -4 1]; + DD_4=(diag(ones(m-2,1),2)-4*diag(ones(m-1,1),1)+6*diag(ones(m,1),0)-4*diag(ones(m-1,1),-1)+diag(ones(m-2,1),-2)); + DD_4(1:2,1:5)=[d4;d4];DD_4(m-1:m,m-4:m)=[d4;d4]; + D = DD_4'*DD_4; + case 5 + d5=0*[-1 5 -10 10 -5 1]; + DD_5=(-diag(ones(m-3,1),-3)+5*diag(ones(m-2,1),-2)-10*diag(ones(m-1,1),-1)+10*diag(ones(m,1),0)-5*diag(ones(m-1,1),1)+diag(ones(m-2,1),2)); + DD_5(1:3,1:6)=[d5;d5;d5]; + DD_5(m-1:m,m-5:m)=[d5;d5]; + D = DD_5'*DD_5; + case 6 + default_arg('scaling', 1/60); + d6=0*[1 -6 15 -20 15 -6 1]; + DD_6=(diag(ones(m-3,1),3)-6*diag(ones(m-2,1),2)+15*diag(ones(m-1,1),1)-20*diag(ones(m,1),0)+15*diag(ones(m-1,1),-1)-6*diag(ones(m-2,1),-2)+diag(ones(m-3,1),-3)); + DD_6(1:3,1:7)=[d6;d6;d6];DD_6(m-2:m,m-6:m)=[d6;d6;d6]; + D = DD_6'*DD_6; + case 7 + d7=0*[-1 7 -21 35 -35 21 -7 1]; + DD_7=(-diag(ones(m-4,1),-4)+7*diag(ones(m-3,1),-3)-21*diag(ones(m-2,1),-2)+35*diag(ones(m-1,1),-1)-35*diag(ones(m,1),0)+21*diag(ones(m-1,1),1)-7*diag(ones(m-2,1),2)+diag(ones(m-3,1),3)); + DD_7(1:4,1:8)=[d7;d7;d7;d7]; + DD_7(m-2:m,m-7:m)=[d7;d7;d7]; + D = DD_7'*DD_7; + case 9 + d9=0*[-1 9 -36 84 -126 126 -84 36 -9 1]; + DD_9=(-diag(ones(m-5,1),-5)+9*diag(ones(m-4,1),-4)-36*diag(ones(m-3,1),-3)+84*diag(ones(m-2,1),-2)-126*diag(ones(m-1,1),-1)+126*diag(ones(m,1),0)-84*diag(ones(m-1,1),1)+36*diag(ones(m-2,1),2)-9*diag(ones(m-3,1),3)+diag(ones(m-4,1),4)); + DD_9(1:5,1:10)=[d9;d9;d9;d9;d9]; + DD_9(m-3:m,m-9:m)=[d9;d9;d9;d9]; + D = DD_9'*DD_9; + otherwise + error('Order not yet supported', order); + end + scaling + D = scaling*Hinv*D; +end
--- a/+scheme/Burgers1D.m Tue Sep 11 16:29:00 2018 +0200 +++ b/+scheme/Burgers1D.m Thu Sep 13 14:21:47 2018 +0200 @@ -6,7 +6,6 @@ params D % Non-stabalized scheme operator - M % Derivative norm H % Discrete norm Hi % Norm inverse e_l @@ -16,44 +15,63 @@ end methods - function obj = Burgers1D(grid, pde_form, operator_type, order, params) + function obj = Burgers1D(grid, pde_form, operator_type, order, dissipation, params) assert(grid.D == 1); assert(grid.size() == length(params.eps)); m = grid.size(); - h = grid.scaling(); lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids. default_arg('pde_form','skew-symmetric'); default_arg('operator_type','narrow'); + default_arg('dissipation','on'); switch operator_type case 'narrow' ops = sbp.D4Variable(m, lim, order); D1 = ops.D1; - D2 = ops.D2; + D2 = ops.D2; + DissipationOp = -1*sbp.dissipationOperator(m, order, ops.HI); + d_l = ops.d1_l'; + d_r = ops.d1_r'; + case 'upwind' + ops = sbp.D1Upwind(m, lim, order); + D1 = (ops.Dp + ops.Dm)/2; + %D2eps = @(eps) ops.Dp*diag(eps)*ops.Dm; + %D2eps = @(eps) ops.Dm*diag(eps)*ops.Dp; + D2 = @(eps) (ops.Dp*diag(eps)*ops.Dm + ops.Dm*diag(eps)*ops.Dp)/2; + DissipationOp = (ops.Dp-ops.Dm)/2; + d_l = D1; + d_r = D1; otherwise error('Other operator types not yet supported', operator_type); end + %% TODO: Figure out how to evaluate viscosity as viscosity(v,t) here instead of parametrizing D on the viscosity. switch pde_form case 'skew-symmetric' - D = @(v, viscosity) -1/3*v.*D1*v - 1/3*D1*v.^2 + D2(params.eps + viscosity)*v; + if (strcmp(dissipation,'on')) + D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v; + else + D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity))*v; + end case 'conservative' - D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v; + if (strcmp(dissipation,'on')) + D = @(v, viscosity) -1/2*D1*v.^2 + (D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v; + else + D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v; + end end obj.grid = grid; obj.order = order; obj.params = params; - %% TODO: Figure out how to evaluate viscosity as viscosity(v,t) here instead of parametrizing D on the viscosity. obj.D = D; - obj.M = ops.M; obj.H = ops.H; obj.Hi = ops.HI; obj.e_l = ops.e_l; obj.e_r = ops.e_r; - obj.d_l = ops.d1_l; - obj.d_r = ops.d1_r; + obj.d_l = d_l; + obj.d_r = d_r; end % Closure functions return the opertors applied to the own doamin to close the boundary @@ -71,7 +89,7 @@ % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary case {'R','robin'} p = s*obj.Hi*e; - closure = @(v, viscosity) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*(obj.params.eps + viscosity)*d'*v); + closure = @(v, viscosity) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*((obj.params.eps + viscosity).*d*v)); switch class(data) case 'double' penalty = s*p*data;