Mercurial > repos > public > sbplib

view +sbp/+implementations/d2_8.m @ 1289:2fd2e2337b77 feature/boundary_optimized_grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add utility function for constructing a (possibly multidimensional) grid based on the grid points used by the boundary optimized SBP operators
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 01 Jul 2020 15:15:30 +0200
parents f7ac3cd6eeaa
children
line wrap: on
line source

function [H, HI, D1, D2, e_1, e_m, M, Q, S_1, S_m] = d2_8(m,h)
    

    BP = 8;
    if(m<2*BP)
        error(['Operator requires at least ' num2str(2*BP) ' grid points']);
    end

    H=speye(m,m);
    H(1:8,1:8)=diag([1498139/5080320, 1107307/725760, 20761/80640, 1304999/725760, 299527/725760, 103097/80640, 670091/725760, 5127739/5080320]);
    H(m-7:m,m-7:m)=rot90(diag([1498139/5080320, 1107307/725760, 20761/80640, 1304999/725760, 299527/725760, 103097/80640, 670091/725760, 5127739/5080320]),2);

%     D1=-(1/280*diag(ones(m-4,1),4)-4/105*diag(ones(m-3,1),3)+1/5*diag(ones(m-2,1),2)-4/5*diag(ones(m-1,1),1)+4/5*diag(ones(m-1,1),-1)-1/5*diag(ones(m-2,1),-2)+4/105*diag(ones(m-3,1),-3)-1/280*diag(ones(m-4,1),-4));

    diags   = -4:4;
    stencil = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280];
    D1 = stripeMatrix(stencil, diags, m);


    %r68 = -1022551/30481920;
    %r78 = 6445687/8709120;
    %r67 = 1714837/4354560;

    %r67=0.58;
    %r68=-0.08;
    %r78=0.75;

    %r67=0.65;
    %r68=-0.1;
    %r78=0.75;

    %r67=0.9250;
    %r68=-0.2;
    %r78=0.775;

    %r67=0.65;
    %r68=-0.105;
    %r78=0.755;

    %%r67=0.649;
    %%r68=-0.104;
    %%r78=0.755;

    %r67=-0.48;
    %r68=0.3;
    %r78=0.67;

    %r67=0.5600;
    %r68=-0.0733;
    %r78=0.7500;/scr0/home/ken/VERY_FINE

    % min med 1/10 f?r D6
    %r67=0.62;
    %r68=-0.1040;
    %r78=0.7640;

    % Den nya optimerade, for att fungera i NS-dissipation

    r67=0.69789473684211;
    r68=-0.12052631578947;
    r78=0.75868421052632;


    D1(1:8,1:12)=[-2540160/1498139, -142642467/5992556+50803200/1498139*r78+5080320/1498139*r67+25401600/1498139*r68, 705710031/5992556-228614400/1498139*r78-25401600/1498139*r67-121927680/1498139*r68, -3577778591/17977668+381024000/1498139*r78+50803200/1498139*r67+228614400/1498139*r68, 203718909/1498139-254016000/1498139*r78-50803200/1498139*r67-203212800/1498139*r68, -32111205/5992556+25401600/1498139*r67+76204800/1498139*r68, -652789417/17977668+76204800/1498139*r78-5080320/1498139*r67, 74517981/5992556-25401600/1498139*r78-5080320/1498139*r68, 0, 0, 0, 0;142642467/31004596-7257600/1107307*r78-725760/1107307*r67-3628800/1107307*r68, 0, -141502371/2214614+91445760/1107307*r78+10886400/1107307*r67+50803200/1107307*r68, 159673719/1107307-203212800/1107307*r78-29030400/1107307*r67-127008000/1107307*r68, -1477714693/13287684+152409600/1107307*r78+32659200/1107307*r67+127008000/1107307*r68, 11652351/2214614-17418240/1107307*r67-50803200/1107307*r68, 36069450/1107307-50803200/1107307*r78+3628800/1107307*r67, -536324953/46506894+17418240/1107307*r78+3628800/1107307*r68, 0, 0, 0, 0;-18095129/134148+3628800/20761*r78+403200/20761*r67+1935360/20761*r68, 47167457/124566-10160640/20761*r78-1209600/20761*r67-5644800/20761*r68, 0, -120219461/124566+25401600/20761*r78+4032000/20761*r67+16934400/20761*r68, 249289259/249132-25401600/20761*r78-6048000/20761*r67-22579200/20761*r68, -2611503/41522+3628800/20761*r67+10160640/20761*r68, -7149666/20761+10160640/20761*r78-806400/20761*r67, 37199165/290654-3628800/20761*r78-806400/20761*r68, 0, 0, 0, 0;3577778591/109619916-54432000/1304999*r78-7257600/1304999*r67-32659200/1304999*r68, -159673719/1304999+203212800/1304999*r78+29030400/1304999*r67+127008000/1304999*r68, 360658383/2609998-228614400/1304999*r78-36288000/1304999*r67-152409600/1304999*r68, 0, -424854441/5219996+127008000/1304999*r78+36288000/1304999*r67+127008000/1304999*r68, 22885113/2609998-29030400/1304999*r67-76204800/1304999*r68, 158096578/3914997-76204800/1304999*r78+7257600/1304999*r67, -296462325/18269986+29030400/1304999*r78+7257600/1304999*r68, 0, 0, 0, 0;-203718909/2096689+36288000/299527*r78+7257600/299527*r67+29030400/299527*r68, 1477714693/3594324-152409600/299527*r78-32659200/299527*r67-127008000/299527*r68, -747867777/1198108+228614400/299527*r78+54432000/299527*r67+203212800/299527*r68, 424854441/1198108-127008000/299527*r78-36288000/299527*r67-127008000/299527*r68, 0, -17380335/1198108+10886400/299527*r67+25401600/299527*r68, -67080435/1198108+25401600/299527*r78-3628800/299527*r67, 657798011/25160268-10886400/299527*r78-3628800/299527*r68, -2592/299527, 0, 0, 0;1529105/1237164-403200/103097*r67-1209600/103097*r68, -3884117/618582+1935360/103097*r67+5644800/103097*r68, 2611503/206194-3628800/103097*r67-10160640/103097*r68, -7628371/618582+3225600/103097*r67+8467200/103097*r68, 5793445/1237164-1209600/103097*r67-2822400/103097*r68, 0, 80640/103097*r67, 80640/103097*r68, 3072/103097, -288/103097, 0, 0;93255631/8041092-10886400/670091*r78+725760/670091*r67, -36069450/670091+50803200/670091*r78-3628800/670091*r67, 64346994/670091-91445760/670091*r78+7257600/670091*r67, -158096578/2010273+76204800/670091*r78-7257600/670091*r67, 67080435/2680364-25401600/670091*r78+3628800/670091*r67, -725760/670091*r67, 0, 725760/670091*r78, -145152/670091, 27648/670091, -2592/670091, 0;-3921999/1079524+25401600/5127739*r78+5080320/5127739*r68, 536324953/30766434-121927680/5127739*r78-25401600/5127739*r68, -334792485/10255478+228614400/5127739*r78+50803200/5127739*r68, 296462325/10255478-203212800/5127739*r78-50803200/5127739*r68, -657798011/61532868+76204800/5127739*r78+25401600/5127739*r68, -5080320/5127739*r68, -5080320/5127739*r78, 0, 4064256/5127739, -1016064/5127739, 193536/5127739, -18144/5127739];


    D1(m-7:m,m-11:m)=rot90( -D1(1:8,1:12),2);


    D1=D1/h;

    %DD=-(1/280*diag(ones(m-4,1),4)-4/105*diag(ones(m-3,1),3)+1/5*diag(ones(m-2,1),2)-4/5*diag(ones(m-1,1),1)+4/5*diag(ones(m-1,1),-1)-1/5*diag(ones(m-2,1),-2)+4/105*diag(ones(m-3,1),-3)-1/280*diag(ones(m-4,1),-4));
    %DD(1:4,1:9)=1/280*[-761,2240,-3920,15680/3,-4900,3136,-3920/3,320,-35,
    %                   -35,-446,980,-980,2450/3,-490,196,-140/3,5,
    %                   5,-80,-266,560,-350,560/3,-70,16,-5/3,
    %                   -5/3,20,-140,-126,350,-140,140/3,-10,1];
%     D(m-7:m,m-11:m)=rot90( -D1(1:8,1:12),2);

%     D2=(-1/560*diag(ones(m-4,1),4)+8/315*diag(ones(m-3,1),3)-1/5*diag(ones(m-2,1),2)+8/5*diag(ones(m-1,1),1)+8/5*diag(ones(m-1,1),-1)-1/5*diag(ones(m-2,1),-2)+8/315*diag(ones(m-3,1),-3)-1/560*diag(ones(m-4,1),-4)-205/72*diag(ones(m,1),0));

    diags   = -4:4;
    left_stencil = [-1/560,8/315,-1/5,8/5];
    stencil = [left_stencil,-205/72,fliplr(left_stencil)];
    D2 = stripeMatrix(stencil, diags, m);

    
    D2(1:8,1:12)=[4870382994799/1358976868290, -893640087518/75498714905,926594825119/60398971924, -1315109406200/135897686829,39126983272/15099742981, 12344491342/75498714905, -451560522577/2717953736580, 0, 0, 0, 0, 0;333806012194/390619153855, -154646272029/111605472530, 1168338040/33481641759, 82699112501/133926567036, -171562838/11160547253, -28244698346/167408208795, 11904122576/167408208795, -2598164715/312495323084, 0, 0, 0, 0;7838984095/52731029988, 1168338040/5649753213, -88747895/144865467, 423587231/627750357, -43205598281/22599012852, 4876378562/1883251071, -5124426509/3766502142, 10496900965/39548272491, 0, 0, 0, 0;-94978241528/828644350023, 82699112501/157837019052, 1270761693/13153084921, -167389605005/118377764289, 48242560214/39459254763, -31673996013/52612339684, 43556319241/118377764289, -44430275135/552429566682, 0, 0, 0, 0;1455067816/21132528431, -171562838/3018932633, -43205598281/36227191596, 48242560214/9056797899, -52276055645/6037865266, 57521587238/9056797899, -80321706377/36227191596, 8078087158/21132528431, -1296/299527, 0, 0, 0;10881504334/327321118845, -28244698346/140280479505, 4876378562/9352031967, -10557998671/12469375956, 57521587238/28056095901, -278531401019/93520319670, 73790130002/46760159835, -137529995233/785570685228, 2048/103097, -144/103097, 0, 0;-135555328849/8509847458140, 11904122576/101307707835, -5124426509/13507694378, 43556319241/60784624701, -80321706377/81046166268, 73790130002/33769235945, -950494905688/303923123505, 239073018673/141830790969, -145152/670091, 18432/670091, -1296/670091, 0;0, -2598164715/206729925524, 10496900965/155047444143, -44430275135/310094888286, 425162482/2720130599, -137529995233/620189776572, 239073018673/155047444143, -144648000000/51682481381, 8128512/5127739, -1016064/5127739, 129024/5127739, -9072/5127739];




    D2(m-7:m,m-11:m)=rot90( D2(1:8,1:12) ,2 );

    D2=D2/h^2;

    DS=sparse(m,m);
    DS(1,1:7)=-[-4723/2100, 839/175, -157/35, 278/105, -103/140, -1/175, 6/175];

    DS(m,m-6:m)=fliplr(-[-4723/2100, 839/175, -157/35, 278/105, -103/140, -1/175, 6/175]);
    DS=DS/h;

    H=h*H;
    HI=inv(H);

    %r1=D1*u-u_x;sqrt(r1'*r1)/m
    %r2=D2*u-u_xx;sqrt(r2'*r2)/m

    %te=eig(D1);
    %tm=max(abs(te));
    %plot(real(te),imag(te),'*');
    %grid;
    %xlabel('Real part');
    %ylabel('Imaginary part');
    %title('Spectrum, minimal spectral radius');
    e_1 = sparse(m,1);
    e_1(1)= 1;
    e_m = sparse(m,1);
    e_m(end)= 1;
    S_1 = -DS(1,:)';
    S_m =  DS(end,:)';

    Q = H*D1-(-(e_1*e_1') + (e_m*e_m'));
    M = -(H*D2-(-e_1*S_1' + e_m*S_m'));
end