Mercurial > repos > public > sbplib
changeset 362:ded4156e53e2
Added 2nd order accurate 2nd derivative with variable coefficents in a separate implementation file, used by the class D2Variable.
| author | Martin Almquist <martin.almquist@it.uu.se> |
|---|---|
| date | Tue, 20 Dec 2016 15:04:02 +0100 |
| parents | f18142c1530b |
| children | 8c8fd649bae3 |
| files | +sbp/+implementations/d2_variable_2.m +sbp/D2Variable.m |
| diffstat | 2 files changed, 66 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
diff -r f18142c1530b -r ded4156e53e2 +sbp/+implementations/d2_variable_2.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_variable_2.m Tue Dec 20 15:04:02 2016 +0100 @@ -0,0 +1,59 @@ +function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_2(m,h) + + BP = 1; + if(m<2*BP) + error(['Operator requires at least ' num2str(2*BP) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv(1) = 1/2; + Hv(m:m) = 1/2; + Hv = h*Hv; + H = spdiag(Hv, 0); + HI = spdiag(1./Hv, 0); + + + % Boundary operators + e_l = sparse(m,1); + e_l(1) = 1; + e_r = rot90(e_l, 2); + + d1_l = sparse(m,1); + d1_l(1:3) = 1/h*[-3/2 2 -1/2]; + d1_r = -rot90(d1_l, 2); + + % D1 operator + diags = -1:1; + stencil = [-1/2 0 1/2]; + D1 = stripeMatrix(stencil, diags, m); + + D1(1,1)=-1;D1(1,2)=1;D1(m,m-1)=-1;D1(m,m)=1; + D1(m,m-1)=-1;D1(m,m)=1; + D1=D1/h; + %Q=H*D1 + 1/2*(e_1*e_1') - 1/2*(e_m*e_m'); + + + M=sparse(m,m); + + scheme_width = 3; + scheme_radius = (scheme_width-1)/2; + r = (1+scheme_radius):(m-scheme_radius); + + function D2 = D2_fun(c) + + Mm1 = -c(r-1)/2 - c(r)/2; + M0 = c(r-1)/2 + c(r) + c(r+1)/2; + Mp1 = -c(r)/2 - c(r+1)/2; + + M(r,:) = spdiags([Mm1 M0 Mp1],0:2*scheme_radius,length(r),m); + + + M(1:2,1:2)=[c(1)/2 + c(2)/2 -c(1)/2 - c(2)/2; -c(1)/2 - c(2)/2 c(1)/2 + c(2) + c(3)/2;]; + M(m-1:m,m-1:m)=[c(m-2)/2 + c(m-1) + c(m)/2 -c(m-1)/2 - c(m)/2; -c(m-1)/2 - c(m)/2 c(m-1)/2 + c(m)/2;]; + M=M/h; + + D2=HI*(-M-c(1)*e_l*d1_l'+c(m)*e_r*d1_r'); + end + D2 = @D2_fun; +end \ No newline at end of file
diff -r f18142c1530b -r ded4156e53e2 +sbp/D2Variable.m --- a/+sbp/D2Variable.m Thu Oct 13 09:34:42 2016 +0200 +++ b/+sbp/D2Variable.m Tue Dec 20 15:04:02 2016 +0100 @@ -31,6 +31,13 @@ obj.e_r, obj.d1_l, obj.d1_r] = ... sbp.implementations.d2_variable_4(m,obj.h); obj.borrowing.M.S = 0.2505765857; + case 2 + [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... + obj.e_r, obj.d1_l, obj.d1_r] = ... + sbp.implementations.d2_variable_2(m,obj.h); + obj.borrowing.M.S = 0.3636363636; + % Borrowing const taken from Virta 2014 + otherwise error('Invalid operator order %d.',order); end