2% ALPHARBCPROJ apply projector to Fourier coeff rep for 1D free space BCs
3%
4% [an,bn] = alphaRBCproj(an,bn,phi,phip) applies RBCs to the alpha_n and beta_n
5% vectors (set of Fourier coefficients of u and v=u_t, respectively).
6% phi is a length-w vector giving a (roll-up) window evaluated on
7% the Nyquist grid, and phip the evaluation of its derivative (scaled as if
8% the grid were unit sized). w sets the tolerance, and phi and phip could
9% be as in the history split. The cost is O(N log N), where
10% N=length(an)=length(bn) is the number of coeffs (N/2 the max frequency).
11%
12% Usage:
13% The user's computational domain is Omega = (-pi+2wh,pi-2wh), where
14% h = 2pi/N is the Nyquist grid spacing implied by the number of coefficients,
15% for the periodic domain [-pi,pi). The user simulation must live in Omega.
16% [-pi,-pi+2wh] and [pi-2wh,pi) are padding
17% regions needed to apply the RBCs, and waves in them will be destroyed.
18% Thus the user sacrifices 4w of the N points to the RBCs, in order that
19% the Fourier series well approximates the free-space solution in Omega.
20% To maintain outgoing Fourier coefficients, this RBC projection must be
21% reapplied in tRBC time units, for wh <= tRBC <= 2wh.
22% The lower limit is needed so that traveling waves have cleared the
23% windowing region (otherwise the window hits them multiple times and
24% aliasing error occurs). The upper limit is so that waves do not wrap around
25% into the other transition region, since the assumption is made that the
26% waves in transition region (eg (pi-2wh,pi-wh)) are out-going.
27% Even without this assumption, spectral differentiation would be needed
28% to recover v in the padding regions, but waves would be hit >1 times by
29% the window, causing aliasing error. (For t > 3wh,
30% plain old periodic pollution would also occur back in Omega).
31% The phi and phip should be as in the history split (as given by window.jl).
32%
33% [an,bn,info] = alphaRBCproj(an,bn,phi,phip) also prints & returns debug info.
34%
35% When called with no arguments, a self-test is done.
36%
37% The method is:
38% 1) evaluate u and v on the h grid, via the FFT from the coeffs an and bn.
39% 2) multiplication of u by phi-windows in the padding regions (which
40% broadens the spectral bandwidth, hence gam is needed),
41% 3) building v=u_t=+-u_x on the grid points in the padded regions to give
42% outgoing 1D wave equation solutions there.
43% 4) inverse FFT to the output coefficients an, bn.
44%
45% Notes:
46% i) N even for now.
47% ii) The total padding region could easily be shrunk from 4w to 3w points.
48% The former gives a bit more flexibility in the projection interval.
49% iii) check Fourier series defn, sign, prefactor, may not matter.
51% Barnett 9/8/24
52if nargin==0, test_alphaRBCproj; return; end
53verb = (nargout>2); % verbosity
54N = length(an); assert(length(bn)==N)
55w = length(phi); assert(length(phip)==w)
56assert(4*w<N)
57h = 2*pi/N;
59% 1)
60u0 = fft(an(:)); v0 = fft(bn(:)); % x on [0,2pi), so padding is in middle
62% 2)
63j = N/2-2*w+(1:4*w)'; % indices to bleach out (combine R and L padding)
64phi=phi(:); phip=phip(:); % force col vecs
65phij = [flipud(phi); zeros(2*w,1); phi]; % R then L padding
66pou = ones(N,1); pou(j) = phij; info.pou=pou; % save the POU, full [0,2pi) grid
67u = u0 .* pou;
68v = v0 .* pou; % the naive step 3, alone would cause O(1/w) reflection.
69% This assumes eg, waves in (pi-2wh,pi-wh) are R-going
71% 3) apply prod rule, overwrite v = (phi.u0)_x = phi.u0_x + phi'.u0 in padding,
72% noting that phi.u0_x is already accounted for in v.
73phipj = (1/h) * [flipud(phip); zeros(2*w,1); phip]; % takes d/dx or -d/dx
74v(j) = v(j) + phipj.*u0(j); % add term +- phi'.u0 in prod rule
76% 4)
77an = reshape(ifft(u),size(an)); bn = reshape(ifft(v),size(bn));
80%%%%%% helpers for test
81function [u,v,xg] = show_uv(an,bn) % plot Fourier rep in real space on Nyq grid
82N = length(an);
83xg = -pi + (0:N-1)/N*2*pi;
84u = fftshift(fft(an)); v = fftshift(fft(bn)); % shift since grid starts -pi
85plot(xg, real([u;v]), '.-')
86a=axis; a(1:2)=[-pi,pi]; axis(a); drawnow
88function [an, bn] = propexact(t, kg, an0, bn0) % evol by t the coeffs on k grid
89ak=abs(kg); s=sin(ak*t); c=cos(ak*t); % (taken from tryEFrep.jl)
90sok=s./ak; sok(kg==0) = t; % note all act on vectors
91an = c.*an0 + sok.*bn0;
92bn = -s.*ak.*an0 + c.*bn0;
94%%%%%%
95function test_alphaRBCproj % test can evolve a source-free WE in free space
96addpath ./utils
97verb = 1; % 1,2,.. for debug tests
98tol = 1e-12;
99gam = 0.5; % fractional exceeding of Nyquist of the sigma
101theta = log(1/tol); % set up window on regular grid in Nour's way
102w = ceil(2*theta/(pi*gam));
103phipfun = window(1.0,theta,w); % hack dt=1 (not so happy window needs dt)
104g = (1:w) - 0.5; % not sure of off-by-one shifts here
105phip = phipfun(g)'; % col vec
106phi = blending(phipfun, g, tol);
107if verb>1
108 figure; plot(g, [phi, phip], '-+'); axis tight; xlabel('gridpoints');
109 legend('\phi','\phi'''); drawnow
110 % check that phip is really deriv of phi on the same grid...
111 foldphipe = [phip;flipud(-phip)]; % append flipped so can do periodic diff
112 foldphip = (pi/w) * perispecdiff([phi;flipud(phi)]); % 2pi/2w rescales t
113 fprintf('max err in phip on grid: %.3g\n', norm(foldphip-foldphipe,inf))
114end
116N = 300; assert(mod(N,2)==0) % grid size N even
117x0 = 0.8; % test initial condition center loc
118k0 = N/sqrt(8*theta); % k-width of Gaussian which hits tol by N/2
119kg = [0:N/2-1, -N/2:-1]; % freq k grid, fft ordering (may differ from Nour)
120c0 = 2/(sqrt(2*pi)*k0); % Gaussian height 2 (so splits into packets of 1)
121an = c0 * exp(-0.5*(kg/k0).^2) .* exp(1i*x0*kg); % bump at x0
122bn = 0*an;
123h = 2*pi/N; % implied grid (used inside proj only, plus for plotting)
124tRBC = 1.0*w*h; % may check fails either side of allowed range [1,2] :)
125fprintf('%.3g of N is in Omega; tRBC=%.3g\n',1-4*w/N,tRBC)
126T = 5*pi;
127nt = ceil(T/tRBC)+1; tg=tRBC*(0:nt-1); uxt = nan(N,nt);
128info.pou = nan(1,N);
129if verb, figure; end
130for n=1:nt, t = tg(n); % t step
131 if verb, subplot(2,1,1); [u,v,xg] = show_uv(an,bn);
132 if verb>1, hold on; plot(xg,fftshift(info.pou),'k-'); hold off; end
133 title(sprintf('u and v: t=%6.2f',t));
134 vline([-1 1]*(pi-2*w*h)); vline([-1 1]*(pi-w*h)); % POU on/off pts
135 uxt(:,n) = u(:);
136 subplot(2,1,2); imagesc(xg,tg,log10(abs(uxt))'); axis xy;
137 caxis([-12 0]); xlabel x; ylabel t; colorbar; title('log_{10}|u(x,t)|')
138 vline([-1 1]*(pi-2*w*h)); vline([-1 1]*(pi-w*h)); % POU on/off pts
139 drawnow;
140 end
141 [an,bn,info] = alphaRBCproj(an,bn,phi,phip); % comment to switch off
142 [an, bn] = propexact(tRBC, kg, an, bn); % we plot *after* prop
143 pause(.5); % anim already too slow anyway :(
144end