Mercurial > repos > public > sbplib

diff -r 8f0ee6ba6313 -r 62147d8ce062 +sbp/+implementations/d1_staggered_upwind_2.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d1_staggered_upwind_2.m	Mon Nov 13 11:49:08 2017 -0800
@@ -0,0 +1,73 @@
+function [xp,xm,Hp,Hm,HIp,HIm,Dp,Dm,Dpp,Dpm,Dmp,Dmm] = d1_staggered_upwind_2(n,L)
+
+assert(n-1 >= 6,'Not enough grid points');
+
+np=n;
+nm=n+1;
+h = L/(n-1);
+
+
+%H- norm
+Hm_U=[0.88829e5 / 0.599299e6 0 0; 0 0.8056187e7 / 0.9588784e7 0; 0 0 0.9700117e7 / 0.9588784e7;];
+%H+ norm
+Hp_U=[0.7489557e7 / 0.19177568e8 0 0; 0 0.1386089e7 / 0.1198598e7 0; 0 0 0.18276939e8 / 0.19177568e8;];
+% Upper part of Qpp. Notice that the B matrix is not included here.
+% Qpp+Qpp^T=S/2,Qpp+Qpm^T=0
+Qpp_U=[-0.1e1 / 0.32e2 0.12886609e8 / 0.19177568e8 -0.1349263e7 / 0.9588784e7; -0.10489413e8 / 0.19177568e8 -0.3e1 / 0.16e2 0.16482403e8 / 0.19177568e8; 0.187491e6 / 0.2397196e7 -0.9290815e7 / 0.19177568e8 -0.11e2 / 0.32e2;];
+% Upper part of Qmp. Notice that the B matrix is not included here.
+% Qmp+Qmp^T=S/2,Qmp+Qmm^T=0
+Qmp_U=[-0.2e1 / 0.21e2 0.12495263e8 / 0.16780372e8 -0.7520839e7 / 0.50341116e8; -0.7700871e7 / 0.16780372e8 -0.31e2 / 0.112e3 0.57771939e8 / 0.67121488e8; 0.2726447e7 / 0.50341116e8 -0.31402783e8 / 0.67121488e8 -0.113e3 / 0.336e3;];
+% The staggered + operator, upper part. Notice that the B matrix is not included here.Qp+Qm^T=0
+Qp_U=[-0.801195e6 / 0.2397196e7 0.16507959e8 / 0.19177568e8 -0.509615e6 / 0.19177568e8; -0.219745e6 / 0.1198598e7 -0.2112943e7 / 0.2397196e7 0.2552433e7 / 0.2397196e7; 0.42087e5 / 0.2397196e7 0.395585e6 / 0.19177568e8 -0.19909849e8 / 0.19177568e8;];
+Hp=spdiags(ones(np,1),0,np,np);
+Hp(1:3,1:3)=Hp_U;
+Hp(np-2:np,np-2:np)=fliplr(flipud(Hp_U));
+Hp=Hp*h;
+HIp=inv(Hp);
+
+Hm=spdiags(ones(nm,1),0,nm,nm);
+Hm(1:3,1:3)=Hm_U;
+Hm(nm-2:nm,nm-2:nm)=fliplr(flipud(Hm_U));
+Hm=Hm*h;
+HIm=inv(Hm);
+
+Qpp=spdiags(repmat([-0.3e1 / 0.8e1 -0.3e1 / 0.8e1 0.7e1 / 0.8e1 -0.1e1 / 0.8e1;],[np,1]),-1:2,np,np);
+Qpp(1:3,1:3)=Qpp_U;
+Qpp(np-2:np,np-2:np)=flipud( fliplr(Qpp_U(1:3,1:3) ) )';
+
+Qpm=-Qpp';
+
+Qmp=spdiags(repmat([-0.3e1 / 0.8e1 -0.3e1 / 0.8e1 0.7e1 / 0.8e1 -0.1e1 / 0.8e1;],[nm,1]),-1:2,nm,nm);
+Qmp(1:3,1:3)=Qmp_U;
+Qmp(nm-2:nm,nm-2:nm)=flipud( fliplr(Qmp_U(1:3,1:3) ) )';
+
+Qmm=-Qmp';
+
+
+Bpp=spalloc(np,np,2);Bpp(1,1)=-1;Bpp(np,np)=1;
+Bmp=spalloc(nm,nm,2);Bmp(1,1)=-1;Bmp(nm,nm)=1;
+
+Dpp=HIp*(Qpp+1/2*Bpp) ;
+Dpm=HIp*(Qpm+1/2*Bpp) ;
+
+
+Dmp=HIm*(Qmp+1/2*Bmp) ;
+Dmm=HIm*(Qmm+1/2*Bmp) ;
+
+
+%%% Start with the staggered
+Qp=spdiags(repmat([-1 1],[np,1]),0:1,np,nm);
+Qp(1:3,1:3)=Qp_U;
+Qp(np-2:np,nm-2:nm)=flipud( fliplr(-Qp_U(1:3,1:3) ) );
+Qm=-Qp';
+
+Bp=spalloc(np,nm,2);Bp(1,1)=-1;Bp(np,nm)=1;
+Bm=Bp';
+
+Dp=HIp*(Qp+1/2*Bp) ;
+
+Dm=HIm*(Qm+1/2*Bm) ;
+
+% grids
+xp = h*[0:n]';
+xm = h*[0 1/2+0:n n]';
\ No newline at end of file
diff -r 8f0ee6ba6313 -r 62147d8ce062 +sbp/+implementations/d1_staggered_upwind_4.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d1_staggered_upwind_4.m	Mon Nov 13 11:49:08 2017 -0800
@@ -0,0 +1,76 @@
+function [xp,xm,Hp,Hm,HIp,HIm,Dp,Dm,Dpp,Dpm,Dmp,Dmm] = d1_staggered_upwind_4(n,L)
+
+assert(n-1 >= 8,'Not enough grid points');
+
+np=n;
+nm=n+1;
+h = L/(n-1);
+
+
+%H- norm
+Hm_U=[0.91e2 / 0.720e3 0 0 0; 0 0.325e3 / 0.384e3 0 0; 0 0 0.595e3 / 0.576e3 0; 0 0 0 0.1909e4 / 0.1920e4;];
+%H+ norm
+Hp_U=[0.805e3 / 0.2304e4 0 0 0; 0 0.955e3 / 0.768e3 0 0; 0 0 0.677e3 / 0.768e3 0; 0 0 0 0.2363e4 / 0.2304e4;];
+% Upper part of Qpp. Notice that the B matrix is not included here.
+% Qpp+Qpp^T=S/2,Qpp+Qpm^T=0
+Qpp_U=[-0.11e2 / 0.1536e4 0.1493e4 / 0.2304e4 -0.571e3 / 0.4608e4 -0.13e2 / 0.768e3; -0.697e3 / 0.1152e4 -0.121e3 / 0.1536e4 0.97e2 / 0.128e3 -0.473e3 / 0.4608e4; 0.373e3 / 0.4608e4 -0.139e3 / 0.256e3 -0.319e3 / 0.1536e4 0.1999e4 / 0.2304e4; 0.1e1 / 0.32e2 -0.121e3 / 0.4608e4 -0.277e3 / 0.576e3 -0.143e3 / 0.512e3;];
+
+% Upper part of Qmp. Notice that the B matrix is not included here.
+% Qmp+Qmp^T=S/2,Qmp+Qmm^T=0
+Qmp_U=[-0.209e3 / 0.5970e4 0.13439e5 / 0.17910e5 -0.13831e5 / 0.47760e5 0.10637e5 / 0.143280e6; -0.44351e5 / 0.71640e5 -0.20999e5 / 0.152832e6 0.230347e6 / 0.229248e6 -0.70547e5 / 0.254720e6; 0.3217e4 / 0.15920e5 -0.86513e5 / 0.114624e6 -0.15125e5 / 0.76416e5 0.1087241e7 / 0.1146240e7; -0.1375e4 / 0.28656e5 0.36117e5 / 0.254720e6 -0.655601e6 / 0.1146240e7 -0.211717e6 / 0.764160e6;];
+
+% The staggered + operator, upper part. Notice that the B matrix is not included here.Qp+Qm^T=0
+Qp_U=[-0.338527e6 / 0.1004160e7 0.4197343e7 / 0.4819968e7 0.1423e4 / 0.803328e6 -0.854837e6 / 0.24099840e8; -0.520117e6 / 0.3012480e7 -0.492581e6 / 0.535552e6 0.2476673e7 / 0.2409984e7 0.520117e6 / 0.8033280e7; -0.50999e5 / 0.3012480e7 0.117943e6 / 0.1606656e7 -0.2476673e7 / 0.2409984e7 0.2712193e7 / 0.2677760e7; 0.26819e5 / 0.1004160e7 -0.117943e6 / 0.4819968e7 -0.1423e4 / 0.803328e6 -0.26119411e8 / 0.24099840e8;];
+
+Hp=spdiags(ones(np,1),0,np,np);
+Hp(1:4,1:4)=Hp_U;
+Hp(np-3:np,np-3:np)=fliplr(flipud(Hp_U));
+Hp=Hp*h;
+HIp=inv(Hp);
+
+Hm=spdiags(ones(nm,1),0,nm,nm);
+Hm(1:4,1:4)=Hm_U;
+Hm(nm-3:nm,nm-3:nm)=fliplr(flipud(Hm_U));
+Hm=Hm*h;
+HIm=inv(Hm);
+
+Qpp=spdiags(repmat([0.7e1 / 0.128e3 -0.67e2 / 0.128e3 -0.55e2 / 0.192e3 0.61e2 / 0.64e2 -0.29e2 / 0.128e3 0.11e2 / 0.384e3;],[np,1]),-2:3,np,np);
+Qpp(1:4,1:4)=Qpp_U;
+Qpp(np-3:np,np-3:np)=flipud( fliplr(Qpp_U(1:4,1:4) ) )';
+
+Qpm=-Qpp';
+
+Qmp=spdiags(repmat([0.7e1 / 0.128e3 -0.67e2 / 0.128e3 -0.55e2 / 0.192e3 0.61e2 / 0.64e2 -0.29e2 / 0.128e3 0.11e2 / 0.384e3;],[nm,1]),-2:3,nm,nm);
+Qmp(1:4,1:4)=Qmp_U;
+Qmp(nm-3:nm,nm-3:nm)=flipud( fliplr(Qmp_U(1:4,1:4) ) )';
+
+Qmm=-Qmp';
+
+
+Bpp=spalloc(np,np,2);Bpp(1,1)=-1;Bpp(np,np)=1;
+Bmp=spalloc(nm,nm,2);Bmp(1,1)=-1;Bmp(nm,nm)=1;
+
+Dpp=HIp*(Qpp+1/2*Bpp) ;
+Dpm=HIp*(Qpm+1/2*Bpp) ;
+
+
+Dmp=HIm*(Qmp+1/2*Bmp) ;
+Dmm=HIm*(Qmm+1/2*Bmp) ;
+
+
+%%% Start with the staggered
+Qp=spdiags(repmat([1/24 -9/8 9/8 -1/24],[np,1]),-1:2,np,nm);
+Qp(1:4,1:4)=Qp_U;
+Qp(np-3:np,nm-3:nm)=flipud( fliplr(-Qp_U(1:4,1:4) ) );
+Qm=-Qp';
+
+Bp=spalloc(np,nm,2);Bp(1,1)=-1;Bp(np,nm)=1;
+Bm=Bp';
+
+Dp=HIp*(Qp+1/2*Bp) ;
+
+Dm=HIm*(Qm+1/2*Bm) ;
+
+% grids
+xp = h*[0:n]';
+xm = h*[0 1/2+0:n n]';
\ No newline at end of file
diff -r 8f0ee6ba6313 -r 62147d8ce062 +sbp/+implementations/d1_staggered_upwind_6.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d1_staggered_upwind_6.m	Mon Nov 13 11:49:08 2017 -0800
@@ -0,0 +1,78 @@
+function [xp,xm,Hp,Hm,HIp,HIm,Dp,Dm,Dpp,Dpm,Dmp,Dmm] = d1_staggered_upwind_6(n,L)
+
+assert(n-1 >= 12,'Not enough grid points');
+
+np=n;
+nm=n+1;
+h = L/(n-1);
+
+bp=6; % Number of boundary points
+
+
+%H- norm
+Hm_U=[0.18373e5 / 0.136080e6 0 0 0 0 0; 0 0.228247e6 / 0.276480e6 0 0 0 0; 0 0 0.219557e6 / 0.207360e6 0 0 0; 0 0 0 0.44867e5 / 0.46080e5 0 0; 0 0 0 0 0.487757e6 / 0.483840e6 0; 0 0 0 0 0 0.2485447e7 / 0.2488320e7;];
+%H+ norm
+Hp_U=[0.174529e6 / 0.552960e6 0 0 0 0 0; 0 0.769723e6 / 0.552960e6 0 0 0 0; 0 0 0.172613e6 / 0.276480e6 0 0 0; 0 0 0 0.343867e6 / 0.276480e6 0 0; 0 0 0 0 0.503237e6 / 0.552960e6 0; 0 0 0 0 0 0.560831e6 / 0.552960e6;];
+% Upper part of Qpp. Notice that the B matrix is not included here.
+% Qpp+Qpp^T=S/2,Qpp+Qpm^T=0
+Qpp_U=[-0.11e2 / 0.12288e5 0.3248629e7 / 0.4976640e7 -0.742121e6 / 0.9953280e7 -0.173089e6 / 0.1658880e7 0.100627e6 / 0.9953280e7 0.263e3 / 0.15552e5; -0.3212989e7 / 0.4976640e7 -0.341e3 / 0.20480e5 0.488429e6 / 0.995328e6 0.2108717e7 / 0.9953280e7 0.5623e4 / 0.829440e6 -0.468779e6 / 0.9953280e7; 0.635201e6 / 0.9953280e7 -0.2135641e7 / 0.4976640e7 -0.2233e4 / 0.30720e5 0.298757e6 / 0.497664e6 -0.2148901e7 / 0.9953280e7 0.99581e5 / 0.1658880e7; 0.184969e6 / 0.1658880e7 -0.2671829e7 / 0.9953280e7 -0.1044721e7 / 0.2488320e7 -0.4697e4 / 0.30720e5 0.882599e6 / 0.995328e6 -0.2114063e7 / 0.9953280e7; -0.118447e6 / 0.9953280e7 0.15761e5 / 0.829440e6 0.915757e6 / 0.9953280e7 -0.2923243e7 / 0.4976640e7 -0.4279e4 / 0.20480e5 0.4614511e7 / 0.4976640e7; -0.263e3 / 0.15552e5 0.422447e6 / 0.9953280e7 -0.5185e4 / 0.331776e6 0.424727e6 / 0.9953280e7 -0.570779e6 / 0.995328e6 -0.2761e4 / 0.12288e5;];
+
+% Upper part of Qmp. Notice that the B matrix is not included here.
+% Qmp+Qmp^T=S/2,Qmp+Qmm^T=0
+Qmp_U=[-0.6660404399e10 / 0.975680535935e12 0.42802131970831759e17 / 0.60695134779444480e17 -0.27241603626152813e17 / 0.91042702169166720e17 0.148772145985039e15 / 0.1264481974571760e16 -0.386865438537449e15 / 0.30347567389722240e17 -0.739491571084877e15 / 0.182085404338333440e18; -0.40937739835891039e17 / 0.60695134779444480e17 -0.20934998251893e14 / 0.570912496455680e15 0.3069347264824655e16 / 0.3082927480860672e16 -0.531249064089227e15 / 0.1541463740430336e16 0.6343721417240047e16 / 0.107902461830123520e18 0.194756943144697e15 / 0.138731736638730240e18; 0.24212381292051773e17 / 0.91042702169166720e17 -0.13932096926615587e17 / 0.15414637404303360e17 -0.22149140293797e14 / 0.285456248227840e15 0.443549615179363e15 / 0.481707418884480e15 -0.1531771257444041e16 / 0.5994581212784640e16 0.23579222779798361e17 / 0.416195209916190720e18; -0.119651391006031e15 / 0.1264481974571760e16 0.2061363806050549e16 / 0.7707318702151680e16 -0.355212104634871e15 / 0.481707418884480e15 -0.42909037900311e14 / 0.285456248227840e15 0.25466291778687943e17 / 0.26975615457530880e17 -0.3289076301679109e16 / 0.11560978053227520e17; 0.153994229603129e15 / 0.30347567389722240e17 -0.2655886084488631e16 / 0.107902461830123520e18 0.782300684927837e15 / 0.5994581212784640e16 -0.17511430871269903e17 / 0.26975615457530880e17 -0.413193098349471e15 / 0.1998193737594880e16 0.189367309285289755e18 / 0.194224431294222336e18; 0.894580992211517e15 / 0.182085404338333440e18 -0.209441083772219e15 / 0.27746347327746048e17 -0.4946149632449393e16 / 0.416195209916190720e18 0.334964642443661e15 / 0.2890244513306880e16 -0.604469352802317407e18 / 0.971122156471111680e18 -0.128172128502407e15 / 0.570912496455680e15;];
+
+% The staggered + operator, upper part. Notice that the B matrix is not included here.Qp+Qm^T=0
+Qp_U=[-0.34660470729017653729e20 / 0.113641961250214656000e21 0.351671379135966469961e21 / 0.415604886857927884800e21 0.1819680091728191503e19 / 0.103901221714481971200e21 -0.18252344147469061739e20 / 0.346337405714939904000e21 -0.18145368485798816351e20 / 0.727308552001373798400e21 0.2627410615589536403e19 / 0.138534962285975961600e21; -0.1606450873889019037e19 / 0.7576130750014310400e19 -0.8503979509850519441e19 / 0.9235664152398397440e19 0.2208731907526094393e19 / 0.2308916038099599360e19 0.1143962309827873891e19 / 0.7696386793665331200e19 0.1263616990270014071e19 / 0.16162412266697195520e20 -0.1402288892096389187e19 / 0.27706992457195192320e20; -0.502728075208147729e18 / 0.11364196125021465600e20 0.5831273443201206481e19 / 0.41560488685792788480e20 -0.9031420599281409001e19 / 0.10390122171448197120e20 0.29977986617775158621e20 / 0.34633740571493990400e20 -0.7995649008389734663e19 / 0.72730855200137379840e20 0.728315692435313537e18 / 0.41560488685792788480e20; 0.710308100786581369e18 / 0.11364196125021465600e20 -0.2317346723533341809e19 / 0.41560488685792788480e20 -0.1357359577229545879e19 / 0.10390122171448197120e20 -0.35124499190079631261e20 / 0.34633740571493990400e20 0.81675241511291974823e20 / 0.72730855200137379840e20 0.144034831596315317e18 / 0.13853496228597596160e20; 0.13360631165154733e17 / 0.841792305557145600e18 -0.875389186128426797e18 / 0.27706992457195192320e20 0.95493318392786453e17 / 0.6926748114298798080e19 0.1714625642820850967e19 / 0.23089160380995993600e20 -0.53371483072841696197e20 / 0.48487236800091586560e20 0.30168063964639488547e20 / 0.27706992457195192320e20; -0.1943232250834614071e19 / 0.113641961250214656000e21 0.8999269402554660119e19 / 0.415604886857927884800e21 0.1242786058815312817e19 / 0.103901221714481971200e21 -0.7480218714053301461e19 / 0.346337405714939904000e21 0.28468183824757664191e20 / 0.727308552001373798400e21 -0.476082521721490837529e21 / 0.415604886857927884800e21;];
+
+Hp=spdiags(ones(np,1),0,np,np);
+Hp(1:bp,1:bp)=Hp_U;
+Hp(np-bp+1:np,np-bp+1:np)=fliplr(flipud(Hp_U));
+Hp=Hp*h;
+HIp=inv(Hp);
+
+Hm=spdiags(ones(nm,1),0,nm,nm);
+Hm(1:bp,1:bp)=Hm_U;
+Hm(nm-bp+1:nm,nm-bp+1:nm)=fliplr(flipud(Hm_U));
+Hm=Hm*h;
+HIm=inv(Hm);
+
+Qpp=spdiags(repmat([-0.157e3 / 0.15360e5 0.537e3 / 0.5120e4 -0.3147e4 / 0.5120e4 -0.231e3 / 0.1024e4 0.999e3 / 0.1024e4 -0.1461e4 / 0.5120e4 0.949e3 / 0.15360e5 -0.33e2 / 0.5120e4;],[np,1]),-3:4,np,np);
+Qpp(1:bp,1:bp)=Qpp_U;
+Qpp(np-bp+1:np,np-bp+1:np)=flipud( fliplr(Qpp_U ) )';
+
+Qpm=-Qpp';
+
+Qmp=spdiags(repmat([-0.157e3 / 0.15360e5 0.537e3 / 0.5120e4 -0.3147e4 / 0.5120e4 -0.231e3 / 0.1024e4 0.999e3 / 0.1024e4 -0.1461e4 / 0.5120e4 0.949e3 / 0.15360e5 -0.33e2 / 0.5120e4;],[nm,1]),-3:4,nm,nm);
+Qmp(1:bp,1:bp)=Qmp_U;
+Qmp(nm-bp+1:nm,nm-bp+1:nm)=flipud( fliplr(Qmp_U ) )';
+
+Qmm=-Qmp';
+
+Bpp=spalloc(np,np,2);Bpp(1,1)=-1;Bpp(np,np)=1;
+Bmp=spalloc(nm,nm,2);Bmp(1,1)=-1;Bmp(nm,nm)=1;
+
+
+Dpp=HIp*(Qpp+1/2*Bpp) ;
+Dpm=HIp*(Qpm+1/2*Bpp) ;
+
+
+Dmp=HIm*(Qmp+1/2*Bmp) ;
+Dmm=HIm*(Qmm+1/2*Bmp) ;
+
+
+%%% Start with the staggered
+Qp=spdiags(repmat([-0.3e1 / 0.640e3 0.25e2 / 0.384e3 -0.75e2 / 0.64e2 0.75e2 / 0.64e2 -0.25e2 / 0.384e3 0.3e1 / 0.640e3],[np,1]),-2:3,np,nm);
+Qp(1:bp,1:bp)=Qp_U;
+Qp(np-bp+1:np,nm-bp+1:nm)=flipud( fliplr(-Qp_U ) );
+Qm=-Qp';
+
+Bp=spalloc(np,nm,2);Bp(1,1)=-1;Bp(np,nm)=1;
+Bm=Bp';
+
+Dp=HIp*(Qp+1/2*Bp) ;
+
+Dm=HIm*(Qm+1/2*Bm) ;
+
+% grids
+xp = h*[0:n]';
+xm = h*[0 1/2+0:n n]';
\ No newline at end of file
diff -r 8f0ee6ba6313 -r 62147d8ce062 +sbp/+implementations/d1_staggered_upwind_8.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d1_staggered_upwind_8.m	Mon Nov 13 11:49:08 2017 -0800
@@ -0,0 +1,82 @@
+function [xp,xm,Hp,Hm,HIp,HIm,Dp,Dm,Dpp,Dpm,Dmp,Dmm] = d1_staggered_upwind_8(n, L)
+
+assert(n-1 >= 12,'Not enough grid points');
+
+np=n;
+nm=n+1;
+h = L/(n-1);
+
+bp=8; % Number of boundary points
+
+
+%H- norm
+Hm_U=[0.27058177e8 / 0.194594400e9 0 0 0 0 0 0 0; 0 0.378196601e9 / 0.464486400e9 0 0 0 0 0 0; 0 0 0.250683799e9 / 0.232243200e9 0 0 0 0 0; 0 0 0 0.48820787e8 / 0.51609600e8 0 0 0 0; 0 0 0 0 0.23953873e8 / 0.23224320e8 0 0 0; 0 0 0 0 0 0.55018021e8 / 0.55738368e8 0 0; 0 0 0 0 0 0 0.284772253e9 / 0.283852800e9 0; 0 0 0 0 0 0 0 0.6036097241e10 / 0.6038323200e10;];
+
+%H+ norm
+Hp_U=[0.273925927e9 / 0.928972800e9 0 0 0 0 0 0 0; 0 0.1417489391e10 / 0.928972800e9 0 0 0 0 0 0; 0 0 0.26528603e8 / 0.103219200e9 0 0 0 0 0; 0 0 0 0.66843235e8 / 0.37158912e8 0 0 0 0; 0 0 0 0 0.76542481e8 / 0.185794560e9 0 0 0; 0 0 0 0 0 0.132009637e9 / 0.103219200e9 0 0; 0 0 0 0 0 0 0.857580049e9 / 0.928972800e9 0; 0 0 0 0 0 0 0 0.937663193e9 / 0.928972800e9;];
+
+% Upper part of Qpp. Notice that the B matrix is not included here.
+% Qpp+Qpp^T=S/2,Qpp+Qpm^T=0
+Qpp_U=[-0.5053e4 / 0.43253760e8 0.299126491231e12 / 0.442810368000e12 -0.562745123e9 / 0.8856207360e10 -0.59616938177e11 / 0.398529331200e12 -0.5489464381e10 / 0.177124147200e12 0.6994759151e10 / 0.88562073600e11 0.712609553e9 / 0.181149696000e12 -0.14167e5 / 0.998400e6; -0.74652297589e11 / 0.110702592000e12 -0.995441e6 / 0.302776320e9 0.662285893e9 / 0.2236416000e10 0.2364992837e10 / 0.5535129600e10 0.8893339297e10 / 0.49816166400e11 -0.2166314077e10 / 0.8945664000e10 -0.78813191e8 / 0.2952069120e10 0.88810243397e11 / 0.1992646656000e13; 0.16939159e8 / 0.276756480e9 -0.2302477691e10 / 0.8200192000e10 -0.197067e6 / 0.9175040e7 0.46740527413e11 / 0.88562073600e11 -0.4827794723e10 / 0.11070259200e11 0.155404199e9 / 0.1230028800e10 0.15110006581e11 / 0.295206912000e12 -0.2491856531e10 / 0.88562073600e11; 0.7568509969e10 / 0.49816166400e11 -0.19762404121e11 / 0.44281036800e11 -0.2554553593e10 / 0.5535129600e10 -0.19550057e8 / 0.302776320e9 0.145135201e9 / 0.210862080e9 0.1596117533e10 / 0.11070259200e11 0.104983477e9 / 0.3558297600e10 -0.1087937837e10 / 0.25303449600e11; 0.5282544031e10 / 0.177124147200e12 -0.131784830977e12 / 0.797058662400e12 0.8310607171e10 / 0.22140518400e11 -0.56309573e8 / 0.105431040e9 -0.36477607e8 / 0.302776320e9 0.57443289107e11 / 0.88562073600e11 -0.220576549e9 / 0.691891200e9 0.39942520697e11 / 0.398529331200e12; -0.1743516779e10 / 0.22140518400e11 0.7784403199e10 / 0.32800768000e11 -0.459036521e9 / 0.4920115200e10 -0.5749179391e10 / 0.22140518400e11 -0.562460513e9 / 0.1383782400e10 -0.1500741e7 / 0.9175040e7 0.70986504791e11 / 0.73801728000e11 -0.930160589e9 / 0.3406233600e10; -0.712609553e9 / 0.181149696000e12 0.180761e6 / 0.6589440e7 -0.18016744831e11 / 0.295206912000e12 0.18651112477e11 / 0.797058662400e12 0.7155507361e10 / 0.44281036800e11 -0.49358401541e11 / 0.73801728000e11 -0.55032223e8 / 0.302776320e9 0.38275518139e11 / 0.40255488000e11; 0.14167e5 / 0.998400e6 -0.88810243397e11 / 0.1992646656000e13 0.650307179e9 / 0.22140518400e11 0.474654619e9 / 0.16102195200e11 -0.7267828861e10 / 0.199264665600e12 0.42227833e8 / 0.425779200e9 -0.71245142351e11 / 0.110702592000e12 -0.7998899e7 / 0.43253760e8;];
+
+% Upper part of Qmp. Notice that the B matrix is not included here.
+% Qmp+Qmp^T=S/2,Qmp+Qmm^T=0
+Qmp_U=[-0.92025012754706822244637e23 / 0.73350939131274317328275670e26 0.930302337374620084855601123690977e33 / 0.1359398284732080668181467336976000e34 -0.727851704787797291000113057457371e33 / 0.2718796569464161336362934673952000e34 0.7360201789777281105403766444579e31 / 0.90626552315472044545431155798400e32 0.13716157770579047943179985700303e32 / 0.543759313892832267272586934790400e33 -0.553550686983724599329589830113e30 / 0.21750372555713290690903477391616e32 0.39889728754596585193681401053e29 / 0.20596943708061828305779808136000e32 0.4597826214543803453588826224861e31 / 0.2718796569464161336362934673952000e34; -0.460817194517403953445488602736801e33 / 0.679699142366040334090733668488000e33 -0.2251237342833501172927156567e28 / 0.267062619272621870023659683840e30 0.190175037040902421814031988319053e33 / 0.202800676510147232549216572416000e33 -0.280264041304585424221475951435491e33 / 0.1081603608054118573595821719552000e34 -0.2948540748840639854083293613843e31 / 0.81120270604058893019686628966400e32 0.10373234521893519136055741251159e32 / 0.194688649449741343247247909519360e33 -0.391946441371315598560753809131e30 / 0.59488198442976521547770194575360e32 -0.15471439859462263133402977429037e32 / 0.6026077244872946338605292437504000e34; 0.702728186606917922841148808013121e33 / 0.2718796569464161336362934673952000e34 -0.738964213741742941634289388463587e33 / 0.811202706040588930196866289664000e33 -0.6938440904692701484464706797e28 / 0.267062619272621870023659683840e30 0.13587935221144065448123508546711e32 / 0.16900056375845602712434714368000e32 -0.398581023027193812136121833051e30 / 0.3244810824162355720787465158656e31 -0.4390487184318690815376655433561e31 / 0.486721623624353358118119773798400e33 0.39375020562052974775725982794643e32 / 0.5948819844297652154777019457536000e34 -0.337431206994896862016299739723e30 / 0.1054563517852765609255926176563200e34; -0.1626484714285090813572964936931e31 / 0.22656638078868011136357788949600e32 0.247081446636583523913555235517491e33 / 0.1081603608054118573595821719552000e34 -0.98606689532786656855690277893813e32 / 0.135200451006764821699477714944000e33 -0.54412910204947725598668603049e29 / 0.801187857817865610070979051520e30 0.44527159063104584858762191285733e32 / 0.54080180402705928679791085977600e32 -0.4170461618635408764831166558231e31 / 0.18541776138070604118785515192320e32 0.2733053052508781643445069010989e31 / 0.55081665224978260692379809792000e32 -0.62697761224771731704532896739409e32 / 0.7030423452351770728372841177088000e34; -0.16732447706243527499842078942603e32 / 0.543759313892832267272586934790400e33 0.166114324967929244010272382781e30 / 0.2897152521573531893560236748800e31 0.45101984053799299547057123965e29 / 0.811202706040588930196866289664e30 -0.17968772701532140224971787578029e32 / 0.27040090201352964339895542988800e32 -0.48635228792591105226576208151e29 / 0.400593928908932805035489525760e30 0.89480251714879473157102672858403e32 / 0.97344324724870671623623954759680e32 -0.42143863277048975845235502012773e32 / 0.148720496107441303869425486438400e33 0.4389611982769118812600028148913e31 / 0.52728175892638280462796308828160e32; 0.590475930631333482244268674627e30 / 0.21750372555713290690903477391616e32 -0.1702402353852012839500776925027e31 / 0.27812664207105906178178272788480e32 0.5538920825569170589993365115709e31 / 0.121680405906088339529529943449600e33 0.13866428236712167588433125984777e32 / 0.129792432966494228831498606346240e33 -0.16470253483258635888811378530287e32 / 0.24336081181217667905905988689920e32 -0.391812951622286986994545412081e30 / 0.2403563573453596830212937154560e31 0.345131812155051195787190432571781e33 / 0.356929190657859129286621167452160e33 -0.8126635967231661246920267342728147e34 / 0.25309524428466374622142228237516800e35; -0.88170259036818455465427409231e29 / 0.41193887416123656611559616272000e32 0.943459254925453305858100463659e30 / 0.118976396885953043095540389150720e33 -0.104249970715879362499981962934393e33 / 0.5948819844297652154777019457536000e34 0.78160374875076232344526852587e29 / 0.18360555074992753564126603264000e32 0.37536035579237342566706229296521e32 / 0.297440992214882607738850972876800e33 -0.240812609144054591374890622842541e33 / 0.356929190657859129286621167452160e33 -0.145390363516271219220036634153e30 / 0.801187857817865610070979051520e30 0.38093040928245760105862644889779897e35 / 0.38667328987934739006050626473984000e35; -0.4596580943680048141920105029111e31 / 0.2718796569464161336362934673952000e34 0.15262190991038620305994595494037e32 / 0.6026077244872946338605292437504000e34 0.1782508142749448850000523441843e31 / 0.1054563517852765609255926176563200e34 -0.33510338065544202178903991034341e32 / 0.7030423452351770728372841177088000e34 -0.8226959851063633384855147175189e31 / 0.421825407141106243702370470625280e33 0.3729744974003821244717412110773747e34 / 0.25309524428466374622142228237516800e35 -0.6554525408845613610278033864711443e34 / 0.9666832246983684751512656618496000e34 -0.49383212966905764275823531781e29 / 0.267062619272621870023659683840e30;];
+
+% The staggered + operator, upper part. Notice that the B matrix is not included here.Qp+Qm^T=0
+Qp_U=[-0.6661444046602902130086192779621746771e37 / 0.23342839720855518078613888837079040000e38 0.12367666586033683088530161144944922123649e41 / 0.14797137255429936031548004199961722880000e41 0.4395151478839780503265910147432452357e37 / 0.186518536833150454179176523528929280000e39 -0.87963490365651280757669626389462038003e38 / 0.1345194295948176002868000381814702080000e40 -0.224187853266917333808369723332131873619e39 / 0.5178998039400477611041801469986603008000e40 0.55314264508663245255306547445101748003e38 / 0.2959427451085987206309600839992344576000e40 0.232303691019191510116997121080032335473e39 / 0.7398568627714968015774002099980861440000e40 -0.2816079756494683118054172509551114287e37 / 0.182680706857159704093185237036564480000e39; -0.10353150194859080046224306922028068797e38 / 0.43350988053017390717425793554575360000e38 -0.41810120772068966379630018024826564868789e41 / 0.44391411766289808094644012599885168640000e41 0.14548989645937048285445132077600344947e38 / 0.16118885899161150361163403267932160000e38 0.1353060140543056426043612771323929599939e40 / 0.6341630252327115442092001799983595520000e40 0.52584560489710238297754746193308329991e38 / 0.317081512616355772104600089999179776000e39 -0.322880558646947920788817445027309369183e39 / 0.8878282353257961618928802519977033728000e40 -0.2523050363123545986203815545725359199053e40 / 0.22195705883144904047322006299942584320000e41 0.2171073153815036601208579021528549178587e40 / 0.44391411766289808094644012599885168640000e41; -0.891108828407068615176254357157902263e36 / 0.14450329351005796905808597851525120000e38 0.9207480733237809022786174965361787199767e40 / 0.44391411766289808094644012599885168640000e41 -0.42541158717297365539935073107511004261e38 / 0.62172845611050151393058841176309760000e38 0.3752589973871193536179209598104238527889e40 / 0.4932379085143312010516001399987240960000e40 -0.503978909830206767329311700636686423011e39 / 0.2219570588314490404732200629994258432000e40 -0.25729240076595732139825342259542295873e38 / 0.328825272342887467367733426665816064000e39 0.283553624900223723890865927423934504751e39 / 0.2466189542571656005258000699993620480000e40 -0.129077863120107719239717930927898530811e39 / 0.4035582887844528008604001145444106240000e40; 0.150348887575956035991597139943595613e36 / 0.2000814833216187263881190471749632000e37 -0.20271017027971837857434511014725441049e38 / 0.422775350155141029472800119998906368000e39 -0.28593401740189393456434271848012721991e38 / 0.87041983855470211950282377646833664000e38 -0.3107901953269564253446821657614993689969e40 / 0.2959427451085987206309600839992344576000e40 0.163990240717967340033770840655582198979e39 / 0.147971372554299360315480041999617228800e39 0.1282390462809880153203810439898877805771e40 / 0.5326969411954776971357281511986220236800e40 0.55264142141446425784998699697993554089e38 / 0.1479713725542993603154800419996172288000e40 -0.11464469241412665723941542481946601159e38 / 0.328825272342887467367733426665816064000e39; 0.523088046329055388999216632374482217e36 / 0.8670197610603478143485158710915072000e37 -0.1060078188350823496391913491020744512571e40 / 0.8878282353257961618928802519977033728000e40 0.2489758378222482566046953325470732791e37 / 0.87041983855470211950282377646833664000e38 0.2223322627542429729644405110485179851147e40 / 0.8878282353257961618928802519977033728000e40 -0.405229856535104846314030690968355556777e39 / 0.443914117662898080946440125998851686400e39 0.1490365162783652350291746570529147907183e40 / 0.1775656470651592323785760503995406745600e40 -0.768972334347761129534143069144613163467e39 / 0.4439141176628980809464401259988516864000e40 0.34866718651627637034229483614924933299e38 / 0.1268326050465423088418400359996719104000e40; -0.522920230014765612739540250860177277e36 / 0.14450329351005796905808597851525120000e38 0.1453599225900263767517831305438030179313e40 / 0.44391411766289808094644012599885168640000e41 0.39466073050115777575909309244731209107e38 / 0.435209919277351059751411888234168320000e39 -0.501532616891797921922210650058733239369e39 / 0.4932379085143312010516001399987240960000e40 -0.276688880738967796756224778348832938429e39 / 0.2219570588314490404732200629994258432000e40 -0.346414321340583244146983258114429082007e39 / 0.328825272342887467367733426665816064000e39 0.12835585968364300662881260070481636919e38 / 0.10676145205937904784666669696942080000e38 -0.823464330534329517579044232163467356639e39 / 0.44391411766289808094644012599885168640000e41; -0.52802811745049309885720368328150617e35 / 0.1688999534533145092886719229399040000e37 0.306081878756402097710283710502974342821e39 / 0.4932379085143312010516001399987240960000e40 -0.43281676756908448328630499995327205907e38 / 0.1305629757832053179254235664702504960000e40 -0.136777517026276610492513752852976312597e39 / 0.4932379085143312010516001399987240960000e40 0.7290267185112688983569547959866321207e37 / 0.246618954257165600525800069999362048000e39 0.351639553804386998900032061955116226307e39 / 0.3804978151396269265255201079990157312000e40 -0.2840607921448362125113163437018746783083e40 / 0.2466189542571656005258000699993620480000e40 0.1859197240206073478058984196606793676119e40 / 0.1644126361714437336838667133329080320000e40; 0.5412884950805861225701675068383948563e37 / 0.303456916371121735021980554882027520000e39 -0.1279848124286614098666833963684894093627e40 / 0.44391411766289808094644012599885168640000e41 0.17218770500911534891873213313315117e35 / 0.39564538116122823613764717112197120000e38 0.904771800018341402297740826650985365019e39 / 0.44391411766289808094644012599885168640000e41 0.65377500469171784914136353007658339737e38 / 0.15536994118201432833125404409959809024000e41 -0.211014619053462658139013021424372225649e39 / 0.8878282353257961618928802519977033728000e40 0.195325945159072852492812627815477159003e39 / 0.3170815126163557721046000899991797760000e40 -0.52260858238454846311625894178508086247499e41 / 0.44391411766289808094644012599885168640000e41;];
+
+Hp=spdiags(ones(np,1),0,np,np);
+Hp(1:bp,1:bp)=Hp_U;
+Hp(np-bp+1:np,np-bp+1:np)=fliplr(flipud(Hp_U));
+Hp=Hp*h;
+HIp=inv(Hp);
+
+Hm=spdiags(ones(nm,1),0,nm,nm);
+Hm(1:bp,1:bp)=Hm_U;
+Hm(nm-bp+1:nm,nm-bp+1:nm)=fliplr(flipud(Hm_U));
+Hm=Hm*h;
+HIm=inv(Hm);
+
+tt=[0.1447e4 / 0.688128e6 -0.17119e5 / 0.688128e6 0.8437e4 / 0.57344e5 -0.5543e4 / 0.8192e4 -0.15159e5 / 0.81920e5 0.16139e5 / 0.16384e5 -0.2649e4 / 0.8192e4 0.15649e5 / 0.172032e6 -0.3851e4 / 0.229376e6 0.5053e4 / 0.3440640e7;];
+
+Qpp=spdiags(repmat(tt,[np,1]),-4:5,np,np);
+Qpp(1:bp,1:bp)=Qpp_U;
+Qpp(np-bp+1:np,np-bp+1:np)=flipud( fliplr(Qpp_U ) )';
+
+Qpm=-Qpp';
+
+Qmp=spdiags(repmat(tt,[nm,1]),-4:5,nm,nm);
+Qmp(1:bp,1:bp)=Qmp_U;
+Qmp(nm-bp+1:nm,nm-bp+1:nm)=flipud( fliplr(Qmp_U ) )';
+
+Qmm=-Qmp';
+
+
+Bpp=spalloc(np,np,2);Bpp(1,1)=-1;Bpp(np,np)=1;
+Bmp=spalloc(nm,nm,2);Bmp(1,1)=-1;Bmp(nm,nm)=1;
+
+Dpp=HIp*(Qpp+1/2*Bpp) ;
+Dpm=HIp*(Qpm+1/2*Bpp) ;
+
+
+Dmp=HIm*(Qmp+1/2*Bmp) ;
+Dmm=HIm*(Qmm+1/2*Bmp) ;
+
+
+%%% Start with the staggered
+Qp=spdiags(repmat([0.5e1 / 0.7168e4 -0.49e2 / 0.5120e4 0.245e3 / 0.3072e4 -0.1225e4 / 0.1024e4 0.1225e4 / 0.1024e4 -0.245e3 / 0.3072e4 0.49e2 / 0.5120e4 -0.5e1 / 0.7168e4;],[np,1]),-3:4,np,nm);
+Qp(1:bp,1:bp)=Qp_U;
+Qp(np-bp+1:np,nm-bp+1:nm)=flipud( fliplr(-Qp_U ) );
+Qm=-Qp';
+
+Bp=spalloc(np,nm,2);Bp(1,1)=-1;Bp(np,nm)=1;
+Bm=Bp';
+
+Dp=HIp*(Qp+1/2*Bp) ;
+
+Dm=HIm*(Qm+1/2*Bm) ;
+
+% grids
+xp = h*[0:n]';
+xm = h*[0 1/2+0:n n]';
\ No newline at end of file
diff -r 8f0ee6ba6313 -r 62147d8ce062 +sbp/D1StaggeredUpwind.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/D1StaggeredUpwind.m	Mon Nov 13 11:49:08 2017 -0800
@@ -0,0 +1,95 @@
+classdef D1StaggeredUpwind < sbp.OpSet
+    % Compatible staggered and upwind operators by Ken Mattsson and Ossian O'reilly
+    properties
+        % x_primal: "primal" grid with m points. Equidistant. Called Plus grid in Ossian's paper.
+        % x_dual: "dual" grid with m+1 points. Called Minus grid in Ossian's paper.
+
+        % D1_primal takes FROM dual grid TO primal grid
+        % D1_dual takes FROM primal grid TO dual grid
+
+        D1_primal % SBP operator approximating first derivative
+        D1_dual % SBP operator approximating first derivative
+
+        Dplus_primal  % Upwind operator on primal grid
+        Dminus_primal % Upwind operator on primal grid
+        Dplus_dual % Upwind operator on dual grid
+        Dminus_dual % Upwind operator on dual grid
+
+        H_primal % Norm matrix
+        H_dual % Norm matrix
+        H_primalI % H^-1
+        H_dualI % H^-1
+        e_primal_l % Left boundary operator
+        e_dual_l % Left boundary operator
+        e_primal_r % Right boundary operator
+        e_dual_r % Right boundary operator
+        m % Number of grid points.
+        m_primal % Number of grid points.
+        m_dual % Number of grid points.
+        h % Step size
+        x_primal % grid
+        x_dual % grid
+        x
+        borrowing % Struct with borrowing limits for different norm matrices
+    end
+
+    methods
+        function obj = D1StaggeredUpwind(m,lim,order)
+
+          x_l = lim{1};
+          x_r = lim{2};
+          L = x_r-x_l;
+
+          m_primal = m;
+          m_dual = m+1;
+
+          switch order
+          case 2
+            [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,...
+            obj.H_primalI, obj.H_dualI,...
+            obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,...
+            obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_2(m, L);
+          case 4
+            [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,...
+            obj.H_primalI, obj.H_dualI,...
+            obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,...
+            obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_4(m, L);
+          case 6
+            [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,...
+            obj.H_primalI, obj.H_dualI,...
+            obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,...
+            obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_6(m, L);
+          case 8
+            [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,...
+            obj.H_primalI, obj.H_dualI,...
+            obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,...
+            obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_8(m, L);
+          otherwise
+           error('Invalid operator order %d.',order);
+          end
+
+          obj.m = m;
+          obj.m_primal = m_primal;
+          obj.m_dual = m_dual;
+          obj.h = L/(m-1);
+
+          obj.e_primal_l = sparse(m_primal,1);
+          obj.e_primal_r = sparse(m_primal,1);
+          obj.e_primal_l(1) = 1;
+          obj.e_primal_r(m_primal) = 1;
+
+          obj.e_dual_l = sparse(m_dual,1);
+          obj.e_dual_r = sparse(m_dual,1);
+          obj.e_dual_l(1) = 1;
+          obj.e_dual_r(m_dual) = 1;
+
+          obj.borrowing = [];
+          obj.x = [];
+
+        end
+
+        function str = string(obj)
+            str = [class(obj) '_' num2str(obj.order)];
+        end
+    end
+end