function test_problem(M,N,tc)
% This function is run from the command line taking argumetns
% M: integer - points in discretization in x
% N: interger - points in discretization in y
% tc: 1 or 2 - test case
% 1 - simple product of low order polynomials
% 2 - product of trig functions
% NOTE !!! : You need at LEAST 30 modes in X and Y to
% solve this problem accurately
%
% This m-file solves the 2D Bi-harmonic problem
% [(DX^4+DY4) - alpha*(DX2+DY2) + Beta]*u = f
% with Dirichlet and Neumann boundary conditions using
% 1) Standard Differentiation Matrices and Tau boundary lines
% 2) Quasi-inverse Techniques with Tau lines
% 3) Quasi-inverse Techniques with Galerkin boundary enforcement
%build the 1D grids
x=cheb_grid(M);
y=cheb_grid(N);
%build the 2d grids
[X,Y]=meshgrid(x,y);
% used for building Tau lines
X_VALS=[-1,1];
Y_VALS=[-1,1];
%parameters for the generalized Biharmonic problem
alpha=1; Beta=1;
% set the forcing and analytic solution based on the test case tc
switch tc
case 1
%simple product phi(2,x)*phi(2,y)
u=(256*(-1 + X.^2).^2.*(-1 + 6*X.^2).*(-1 + Y.^2).^2.*(-1 + 6*Y.^2))/25;
f=(256*(Beta*(-1 + X.^2).^2.*(-1 + 6*X.^2).*(-1 + Y.^2).^2.* ...
(-1 + 6*Y.^2) + 24.*(-13 + 90*X.^2).*(-1 + Y.^2).^2.* ...
(-1 + 6*Y.^2) + 24*(-1 + X.^2).^2.*(-1 + 6*X.^2).* ...
(-13 + 90*Y.^2) - 4*alpha*(-8 + 71*Y.^2 - 97*Y.^4 + ...
24*Y.^6 + 6*X.^6.*(4 - 39*Y.^2 + 45*Y.^4) + ...
X.^2.*(71 - 624*Y.^2 + 867*Y.^4 - 234*Y.^6) + ...
X.^4.*(-97 + 867*Y.^2 - 1170*Y.^4 + 270*Y.^6))))/25;
case 2
% u = sin(2*pi*x)*phi_2(x)*sin(2*pi*y)*phi_2(y)
u=(256*(-1 + X.^2).^2.*(-1 + 6*X.^2).*(-1 + Y.^2).^2.*(-1 + 6.*Y.^2).* ...
sin(2*pi*X).*sin(2*pi*Y))/25;
f=(256*(-16*pi*(-1 + X.^2).^2.*(-1 + 6*X.^2).*Y.*...
(4*(39 + alpha + 8*pi^2) - (360 + 13*alpha + 104*pi^2).*Y.^2 + ...
9*(alpha + 8*pi^2).*Y.^4).*cos(2*pi*Y).*sin(2*pi*X) + ...
(-16*pi.*X.*(4*(39 + alpha + 8*pi^2) - ...
(360 + 13*alpha + 104*pi^2).*X.^2 + ...
9*(alpha + 8*pi^2).*X.^4).*(-1 + Y.^2).^2.*(-1 + 6*Y.^2).*...
cos(2*pi*X) + (Beta*(-1 + X.^2).^2.*(-1 + 6*X.^2).* ...
(-1 + Y.^2).^2.*(-1 + 6*Y.^2) + ...
4*alpha*(8 - 71*Y.^2 + 97*Y.^4 - 24*Y.^6 + ...
2*pi^2*(-1 + X.^2).^2.*(-1 + 6*X.^2).*(-1 + Y.^2).^2.* ...
(-1 + 6*Y.^2) - 6*X.^6.*(4 - 39*Y.^2 + 45*Y.^4) + ...
X.^4.*(97 - 867*Y.^2 + 1170*Y.^4 - 270*Y.^6) + ...
X.^2.*(-71 + 624*Y.^2 - 867*Y.^4 + 234*Y.^6)) + ...
8*(4*pi^4.*(-1 + X.^2).^2.*(-1 + 6*X.^2).*(-1 + Y.^2).^2.* ...
(-1 + 6*Y.^2) + 12*pi^2.*(8 - 71*Y.^2 + 97*Y.^4 - ...
24*Y.^6 - 6*X.^6.*(4 - 39*Y.^2 + 45*Y.^4) + ...
X.^4.*(97 - 867*Y.^2 + 1170*Y.^4 - 270*Y.^6) + ...
X.^2.*(-71 + 624*Y.^2 - 867*Y.^4 + 234*Y.^6)) + ...
3*(26 - 194*Y.^2 + 169*Y.^4 - 78*Y.^6 - ...
13*X.^4.*(-13 + 90*Y.^2) + 6*X.^6.*(-13 + 90*Y.^2) + ...
2*X.^2.*(-97 + 45*Y.^2.*(16 - 13*Y.^2 + 6*Y.^4))))).* ...
sin(2*pi*X)).*sin(2*pi*Y)))/25;
end
disp(sprintf('Solving test cases for M = %d and N = %d',M,N))
% Solution type 1 - Full Differentiation Matrices and Tau Boundary Lines
%build the 1D operators
%x operators
%identity matrix
I_X=sparse(eye(M));
%boundary condition identity matrix
I_X_bc=I_X; I_X_bc(end-3:end,:)=0;
%derivative operators
DX4=sparse(deriv(M)^4);
DX2=sparse(deriv(M)^2);
Biharm_X=sparse(DX4-alpha*DX2+.5*Beta*I_X);
%x operators
%identity matrix
I_Y=sparse(eye(N));
%boundary condition identity matrix
I_Y_bc=I_Y; I_Y_bc(end-3:end,:)=0;
%derivative operators
DY4=deriv(N)^4;
DY2=deriv(N)^2;
Biharm_Y=sparse(DY4-alpha*DY2+.5*Beta*I_Y);
% build boundary condition matrice
[u_X_bc,u_X_bc_terms]=bc_2d_new(M,N,'X',{'D_R','D_L','N_R','N_L'},[zeros(N,1),zeros(N,1),zeros(N,1),zeros(N,1)],X_VALS,Y_VALS);
[u_Y_bc,u_Y_bc_terms]=bc_2d_new(M,N,'Y',{'D_R','D_L','N_R','N_L'},[zeros(M,1),zeros(M,1),zeros(M,1),zeros(M,1)],X_VALS,Y_VALS);
%build the "A" Matrix
A=sparse(kron(I_X_bc*Biharm_X,I_Y_bc)+u_X_bc + ... %X operator and X boundary conditions
kron(I_X_bc,I_Y_bc*Biharm_Y)+kron(I_X_bc,I_Y)*u_Y_bc); %Y operator and Y boundary conditions
B=sparse(kron(I_X_bc,I_Y_bc));
figure(1)
spy(A)
title('Full Differential Operators with Tau Lines')
%convert forcing to physical space and stretch
F=phys2cheb2(f); % the 2 at the end of phys2cheb2 signifies a 2 dimensional transform
F=F(:); %stretch
%Solve for U_CALC in spectral space
tic
U_CALC=A\(B*F+u_X_bc_terms+sparse(kron(I_X_bc,I_Y))*u_Y_bc_terms); % since the bcs are zero, we don't really need to add them on here
tts=toc;
%convert back to physical space
U_CALC=reshape(U_CALC,N,M);
u_calc=cheb2phys2(U_CALC);
%determine infinity norm error
full_op_err=max(max(abs(u-u_calc)));
disp(sprintf(' For the Full Differential Operator:'))
disp(sprintf(' Inf-Norm Error : %g',full_op_err))
disp(sprintf(' Time to Solve : %g',tts))
% Solution type 2 - Quasi-Inverse Matrices and Tau Boundary Lines
%build the 1D operators
%x operators
%identity matrix
I_X=sparse(eye(M));
%quasi/pseudo identity
I_4_X=pseudo_eye(M,4);
%quasi-inverse derivative operators
iDX4=sparse(inv_derivP(M,4));
iDX2=sparse(inv_derivP(M,2));
OP_X=sparse(I_4_X*(I_X-alpha*iDX2+.5*Beta*iDX4));
%y operators
%identity matrix
I_Y=sparse(eye(N));
%quasi/pseudo identity matrix
I_4_Y=pseudo_eye(N,4);
%quasi-inverse derivative operators
iDY4=sparse(inv_derivP(N,4));
iDY2=sparse(inv_derivP(N,2));
OP_Y=sparse(I_4_Y*(I_Y-alpha*iDY2+.5*Beta*iDY4));
% build boundary condition matrice
[u_X_bc,u_X_bc_terms]=bc_2d_new_top(M,N,'X',{'D_R','D_L','N_R','N_L'},[zeros(N,1),zeros(N,1),zeros(N,1),zeros(N,1)],X_VALS,Y_VALS);
[u_Y_bc,u_Y_bc_terms]=bc_2d_new_top(M,N,'Y',{'D_R','D_L','N_R','N_L'},[zeros(M,1),zeros(M,1),zeros(M,1),zeros(M,1)],X_VALS,Y_VALS);
%build the "A" Matrix
A=sparse(kron(OP_X,iDY4)+u_X_bc + ... %X operator and X boundary conditions
kron(iDX4,OP_Y)+kron(I_4_X,I_Y)*u_Y_bc); %Y operator and Y boundary conditions
B=sparse(kron(iDX4,iDY4));
figure(2)
subplot(1,2,1)
spy(A)
title('Quasi-Inverse Method LHS Operator A with Tau Lines')
subplot(1,2,2)
spy(B)
title('Quasi-Inverse Method RHS Operator B')
%convert forcing to physical space and stretch
F=phys2cheb2(f); % the 2 at the end of phys2cheb2 signifies a 2 dimensional transform
F=F(:); %stretch
%Solve for U_CALC in spectral space
tic
U_CALC=A\(B*F);
tts=toc;
%convert back to physical space
U_CALC=reshape(U_CALC,N,M);
u_calc=cheb2phys2(U_CALC);
%determine infinity norm error
full_op_err=max(max(abs(u-u_calc)));
disp(sprintf(' For the Quasi-Inverse Technique with Tau Lines:'))
disp(sprintf(' Inf-Norm Error : %g',full_op_err))
disp(sprintf(' Time to Solve : %g',tts))
% Solution type 3 - Quasi-Inverse Matrices and Galerkin Basis Functions
%build the 1D operators
%x operators
%identity matrix
I_X=sparse(eye(M));
%quasi/pseudo identity matrix
I_4_X=pseudo_eye(M,4);
%quasi-inverse derivative operators
iDX4=sparse(inv_derivP(M,4));
iDX2=sparse(inv_derivP(M,2));
OP_X=sparse(I_4_X*(I_X-alpha*iDX2+.5*Beta*iDX4));
%y operators
%identity matrix
I_Y=sparse(eye(N));
%quasi/pseudo identity matrix
I_4_Y=pseudo_eye(N,4);
%quasi-inverse derivative operators
iDY4=sparse(inv_derivP(N,4));
iDY2=sparse(inv_derivP(N,2));
OP_Y=sparse(I_4_Y*(I_Y-alpha*iDY2+.5*Beta*iDY4));
%build stencil matrices
S_X=stencil_mat(M,4);
S_Y=stencil_mat(N,4);
%build the "A" Matrix
A=sparse(kron(restrictA(OP_X*S_X,4),restrictA(iDY4*S_Y,4)) + ...
kron(restrictA(iDX4*S_X,4),restrictA(OP_Y*S_Y,4)));
B=sparse(kron(restrictA(iDX4,4),restrictA(iDY4,4)));
%restrict forcing
F=phys2cheb2(f);
F=F(1:end-4,1:end-4);
F=F(:);
figure(3)
subplot(1,2,1)
spy(A)
title('Restricted Quasi-Inverse Method LHS Operator A with Galerkin Method')
subplot(1,2,2)
spy(B)
title('Restricted Quasi-Inverse Method RHS Operator B')
%Solve for U_CALC in spectral space
tic
U_CALC=A\(B*F);
tts=toc;
%reshape and convert from Galerkin to Chebyshev
U_CALC=reshape(U_CALC,N-4,M-4);
U_CALC=galerkin2cheb2(U_CALC,4);
u_calc=cheb2phys2(U_CALC);
%determine infinity norm error
full_op_err=max(max(abs(u-u_calc)));
disp(sprintf(' For the Quasi-Inverse Technique with Galerkin Method:'))
disp(sprintf(' Inf-Norm Error : %g',full_op_err))
disp(sprintf(' Time to Solve : %g',tts))
