driver_tfe.m

% driver_tfe_small.m
% Script to test the TFE method for DNOs (2d): Small perturbation.

clear all; close all; SavePlots = 0;

% alpha, gamma are the underscore values
alpha = 0; gamma = 1.21; gammaw = 1.15; Eps_Max = 1e-3; Delta_Max = 1e-5;
d = 2*pi; a = 1; Nx = 48; Nz = 32; N = 8; M = 8; b = 1; Test_Case = 17;
Upper_Field = true;
k = sqrt(alpha^2 + gamma^2);
kw = sqrt(alpha^2 + gammaw^2);

% Pass through underscore values in the upper and lower fields and DNOs
x = (d/Nx)*[0:Nx-1]'; p = (2*pi/d)*[0:Nx/2-1,-Nx/2:-1]';
alphap = alpha + p; gammap = 0*alphap; gammapw= 0*alphap;
for j=1:Nx
  if(alphap(j)^2 < k^2)
    gammap(j) = sqrt(k^2 - alphap(j)^2);
  else
    gammap(j) = 1i*sqrt(alphap(j)^2 - k^2);
  end
  
  if(alphap(j)^2 < kw^2)
    gammapw(j) = sqrt(kw^2 - alphap(j)^2);
  else
    gammapw(j) = 1i*sqrt(alphap(j)^2 - kw^2);
  end
end
[Dz,z] = cheb(Nz);
z_orig = z;

% Transform z from [-1,1] to [b=0,a] in upper field
z_min = 0;
z_max = a;
z = ((z_max-z_min)/2.0)*(z - 1.0) + z_max;

% Transform z from [-1,1] to [b=-b,0] in lower field (b is positive constant)
z_min_lf = -b;
z_max_lf = 0;
z_lf = ((z_max_lf-z_min_lf)/2.0)*(z_orig - 1.0) + z_max_lf;

% Keep a copy of the actual values for computing real values with the infinity norm
% Convention is to prefix actual values with "act_"
act_alpha = (1.0 + Delta_Max)*alpha;
act_gamma = (1.0 + Delta_Max)*gamma;
act_gammaw = (1.0 + Delta_Max)*gammaw;
act_k = sqrt(act_alpha^2 + act_gamma^2);
act_kw = sqrt(act_alpha^2 + act_gammaw^2);
act_p = p; act_alphap = act_alpha + p; act_gammap = 0*act_alphap; 
act_gammapw = 0*act_alphap;
for j=1:Nx
  if(act_alphap(j)^2 < act_k^2)
    act_gammap(j) = sqrt(act_k^2 - act_alphap(j)^2);
  else
    act_gammap(j) = 1i*sqrt(act_alphap(j)^2 - act_k^2);
  end
  
  if(act_alphap(j)^2 < act_kw^2)
    act_gammapw(j) = sqrt(act_kw^2 - act_alphap(j)^2);
  else
    act_gammapw(j) = 1i*sqrt(act_alphap(j)^2 - act_kw^2);
  end
end

f = exp(cos(x));
f_x = -sin(x).*f;

% xi and nu should have the actual values
Ar = -3.0; Aw = -4.8; r = 2; alphap_r = alphap(r+1); gammap_r = gammap(r+1);
p_r = p(r+1); gammap_rw = gammapw(r+1); act_alphap_r = act_alphap(r+1); 
act_p_r = act_p(r+1); act_gammap_r = act_gammap(r+1); 
act_gammap_rw = act_gammapw(r+1);

% Make x_full and z_full to pass into w_actual
x_full = zeros(Nx,Nz+1);
f_full = zeros(Nx,Nz+1);
fx_full = zeros(Nx,Nz+1);
z_full = zeros(Nx,Nz+1);
z_full_lf = zeros(Nx,Nz+1);
for ell=0:Nz
  x_full(:,ell+1) = x;
  f_full(:,ell+1) = f;
  fx_full(:,ell+1) = f_x;
end
% Upper Field
for j=1:Nx
  z_full(j,:) = z;
end
% Lower field
for j=1:Nx
  z_full_lf(j,:) = z_lf;
end

% Upper Field
xi = Ar*exp(1i*act_p_r*x).*exp(1i*act_gammap_r*Eps_Max*f);
% nu is the exact solution of the DNO
nu = (-1i*act_gammap_r+1i*act_p_r*Eps_Max*f_x).*xi;
z_pre_full = ((a-Eps_Max*f_full)/a).*z_full + Eps_Max*f_full;
% u is the exact solution of the Upper Field expansion
u = Ar*exp(1i*act_p_r*x_full).*exp(1i*act_gammap_r*z_pre_full);

% Lower Field
xi_lf = Aw*exp(1i*act_p_r*x).*exp(-1i*act_gammap_rw*Eps_Max*f);
% nu is the exact solution of the DNO
nu_lf = (1i*act_gammap_rw+1i*act_p_r*Eps_Max*f_x).*xi_lf;
z_pre_full_lf = ((b+Eps_Max*f_full)/b).*z_full_lf + Eps_Max*f_full;
% w is the exact solution of the Lower Field expansion (multiply gammap by -1)
w = Aw*exp(1i*act_p_r*x_full).*exp(-1i*act_gammap_rw*z_pre_full_lf);

unm = field_tfe_helmholtz_m_and_n(xi,f,p,gammap,alpha,gamma,...
        Dz,a,Nx,Nz,N,M);
wnm = field_tfe_helmholtz_m_and_n_lf(xi_lf,f,p,gammapw,alpha,...
        gammaw,Dz,b,Nx,Nz,N,M);
Gnm_tfe = dno_tfe_helmholtz_m_and_n(unm,f,p,Dz,a,Nx,...
        Nz,N,M);
Gnm_tfe_lf = dno_tfe_helmholtz_m_and_n_lf(wnm,f,p,Dz,...
        b,Nx,Nz,N,M);


% Setup Zeta_nm and Psi_nm for Two Layer Solve (MMS)
[Unm, Wnm] = setup_zeta_and_psi(Ar,Aw,p_r,...
    N,Nx,M,x,f,f_x,xi,xi_lf,nu,nu_lf,Eps_Max,Delta_Max,gammap,...
    gammapw,alpha,gamma,gammaw,Dz,a,Nz,p,r,k,kw,Test_Case,...
    act_gammap,act_gammapw);
                
% Upper Field
if (Upper_Field == true)
  if (N > 0 && M == 0)
    %n_convergence_mms(Unm,Wnm,xi,xi_lf,Eps_Max,N,Nx);
    [relerr,nplot] = n_convergence(u,unm,nu,Gnm_tfe,Eps_Max,N,Nx,Nz);
  elseif (N == 0 && M > 0)
    %m_convergence_mms(Unm,Wnm,xi,xi_lf,Delta_Max,M,Nx);  
    [relerr,nplot] = m_convergence(u,unm,nu,Gnm_tfe,Delta_Max,M,Nx,Nz);
  else
    %nm_convergence_mms(Unm,Wnm,xi,xi_lf,Eps_Max,Delta_Max,N,M,...
    %    Test_Case,Nx);
    %[relerr,nplot,mplot] = nm_convergence(u,unm,nu,Gnm_tfe,Eps_Max,...
    %    Delta_Max,N,M,Nx,Nz);
    make_contour_plots_nm(u,unm,nu,Gnm_tfe,Eps_Max,Delta_Max,N,M,...
        Nx,Nz,Test_Case);
    setup_contour_plots_eps_delta(d,alpha,gamma,Eps_Max,Delta_Max,a,...
        Nx,Nz,N,M,k,Test_Case,r,Ar,f,f_x);
  end
else
% Lower Field
  if (N > 0 && M == 0)
    n_convergence_mms(Unm,Wnm,xi,xi_lf,Eps_Max,N,Nx);
    [relerr,nplot] = n_convergence(w,wnm,nu_lf,Gnm_tfe_lf,Eps_Max,N,Nx,Nz);
  elseif (N == 0 && M > 0)
    m_convergence_mms(Unm,Wnm,xi,xi_lf,Delta_Max,M,Nx);  
    [relerr,nplot] = m_convergence(w,wnm,nu_lf,Gnm_tfe_lf,Delta_Max,...
        M,Nx,Nz);
  else
    nm_convergence_mms(Unm,Wnm,xi,xi_lf,Eps_Max,Delta_Max,N,M,...
        Test_Case,Nx);
    [relerr,nplot,mplot] = nm_convergence(w,wnm,nu_lf,Gnm_tfe_lf,Eps_Max,...
        Delta_Max,N,M,Nx,Nz);
    make_contour_plots_nm(w,wnm,nu_lf,Gnm_tfe_lf,Eps_Max,Delta_Max,...
        N,M,Nx,Nz,Test_Case);
    setup_contour_plots_eps_delta_lf(d,alpha,gammaw,Eps_Max,Delta_Max,b,...
        Nx,Nz,N,M,kw,Test_Case,r,Aw,f,f_x);
  end
end

return;