Mercurial > repos > public > sbplib

view +sbp/+implementations/d1_upwind_4.m @ 774:66eb4a2bbb72 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove default scaling of the system. The scaling doens't seem to help actual solutions. One example that fails in the flexural code. With large timesteps the solutions seems to blow up. One particular example is profilePresentation on the tdb_presentation_figures branch with k = 0.0005
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 18 Jul 2018 15:42:52 -0700
parents 88584b0cfba1
children
line wrap: on
line source

function [H, HI, Dp, Dm, e_1, e_m] = d1_upwind_4(m,h)

    if(m<8)
        error('Operator requires at least 8 grid points');
    end

    Hv=ones(m,1);
    Hv(1:4)=[49/144 61/48 41/48 149/144];
    Hv(m-3:m)=rot90(Hv(1:4),2);
    Hv = Hv*h;
    H = spdiag(Hv,0);
    HI = spdiag(1./Hv,0);

    q_diags   = [-1, 0, 1, 2 3];
    q_stencil = [-1/4 -5/6 3/2 -1/2 1/12];
    Qp = stripeMatrix(q_stencil, q_diags,m);

    Q_U = [
        -0.1e1/0.48e2    0.205e3/0.288e3 -0.29e2/0.144e3 0.1e1/0.96e2;
        -0.169e3/0.288e3 -0.11e2/0.48e2  0.33e2/0.32e2   -0.43e2/0.144e3;
        0.11e2/0.144e3   -0.13e2/0.32e2  -0.29e2/0.48e2  0.389e3/0.288e3;
        0.1e1/0.32e2     -0.11e2/0.144e3 -0.65e2/0.288e3 -0.13e2/0.16e2;
    ];

    Qp(1:4,1:4)=Q_U;
    Qp(m-3:m,m-3:m)=rot90(Q_U, 2)'; %%% This is different from standard SBP

    Qm=-Qp';

    e_1=sparse(m,1);e_1(1)=1;
    e_m=sparse(m,1);e_m(m)=1;

    Dp=HI*(Qp-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;

    Dm=HI*(Qm-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
end