WENO 1D

function out= weno5_DBC(xx,uu,BC1,BC2,BC3,BC4,chose)
% 施加 Diriclet  Boundary  Condition
% 检好！
global x u h
     x=xx;
     u=uu;
     h=x(2)-x(1);
ghostL=[BC1,BC2];
ghostR=[BC3,BC4];
%*****************************************************************
%*********************  1.u^{-}_{i+1/2} **************************
umm = circshift(u,[0 2]);     % v(i-2)
umm(1:2)=ghostL;
um  = circshift(u,[0 1]);     % v(i-1)
um(1)=ghostL(2);
up  = circshift(u,[0 -1]);    % v(i+1)
up(end)=ghostR(1);
upp = circshift(u,[0 -2]);    % v(i+2)
upp(end-1:end)=ghostR;
%2. Polynomials
p0n = (2*umm - 7*um + 11*u)/6;
p1n = ( -um  + 5*u  + 2*up)/6;
p2n = (2*u   + 5*up - upp )/6;
% 3.Smooth Indicators (Beta factors)
B0n = 13/12*(umm-2*um+u  ).^2 + 1/4*(umm-4*um+3*u).^2; 
B1n = 13/12*(um -2*u +up ).^2 + 1/4*(um-up).^2;
B2n = 13/12*(u  -2*up+upp).^2 + 1/4*(3*u-4*up+upp).^2;
%4. Constants
d0n = 1/10; d1n = 6/10; d2n = 3/10; epsilon = 1e-6;
%5. Alpha weights 
alpha0n = d0n./(epsilon + B0n).^2;
alpha1n = d1n./(epsilon + B1n).^2;
alpha2n = d2n./(epsilon + B2n).^2;
alphasumn = alpha0n + alpha1n + alpha2n;
%6. ENO stencils weigths
w0n = alpha0n./alphasumn;
w1n = alpha1n./alphasumn;
w2n = alpha2n./alphasumn;
% Numerical Flux at cell boundary, $u_{i+1/2}^{-}$;
hn = w0n.*p0n + w1n.*p1n + w2n.*p2n;
%*********************************************************************
%************************* u^{+}_{i-1/2} *****************************
  p0p = ( -umm + 5*um + 2*u  )/6;
  p1p = ( 2*um + 5*u  - up   )/6;
  p2p = (11*u  - 7*up + 2*upp)/6;
  B0p = 13/12*(umm-2*um+u  ).^2 + 1/4*(umm-4*um+3*u).^2;  
  B1p = 13/12*(um -2*u +up ).^2 + 1/4*(um-up).^2;   
  B2p = 13/12*(u  -2*up+upp).^2 + 1/4*(3*u -4*up+upp).^2;
% Constants
d0p = 3/10; d1p = 6/10; d2p = 1/10; epsilon = 1e-6;
% Alpha weights 
alpha0p = d0p./(epsilon + B0p).^2;
alpha1p = d1p./(epsilon + B1p).^2;
alpha2p = d2p./(epsilon + B2p).^2;
alphasump = alpha0p + alpha1p + alpha2p;
% ENO stencils weigths
w0p = alpha0p./alphasump;
w1p = alpha1p./alphasump;
w2p = alpha2p./alphasump;
% Numerical Flux at cell boundary, $u_{i-1/2}^{+}$;
hp = w0p.*p0p + w1p.*p1p + w2p.*p2p;
hp2=circshift(hp,[0,-1]);  % u^(+)_{i+1/2};
HP2=hp(end-4:end);
           hp2(end)=La(x(end-4)-h/2,HP2,x(end)+h/2);
hn2=circshift(hn,[0,1]) ; % u^(-)_{i-1/2};
HN2=hn(1:5);
           hn2(1) =La(x(1)+h/2, HN2,x(1)-h/2);
%=============================================================================
% by Prof.Shu Chi-Wang's paper, the left derive of du^{-}=(u^{-}_{i+1/2}-u^{-}_{i-1/2})/dx 
%=============================================================================
if chose==1
    out=(hp2-hp)/h;  % du^+
elseif chose==-1
    out=(hn-hn2)/h;   % du^-
end

    function out=La(x1,y,x)
       out= ((h - x + x1).*(2*h - x + x1)*(3*h - x + x1).*(4*h - x + x1)*y(1) -...
               4*(x - x1).*(-4*h + x - x1).*(-3*h + x - x1).*(-2*h + x - x1)*y(2) +... 
               6*(x - x1).*(-4*h + x - x1).*(-3*h + x - x1).*(-h + x - x1)*y(3) -...
               4*(x - x1).*(-4*h + x - x1).*(-2*h + x - x1).*(-h + x - x1).*y(4) +... 
               (x - x1).*(-3*h + x - x1).*(-2*h + x - x1).*(-h + x - x1).*y(5))./(24.*power(h,4));
    end

end

WENO 1D

标签：lse   ted   als   strong   his   power   sch   size   out   

原文地址：https://www.cnblogs.com/yuewen-chen/p/11587998.html