Create SolveSpectralBIE.m

Solver/SolveSpectralBIE.m

@@ -0,0 +1,43 @@
+function [u, au, bu, cu, FLAG, RELRES, ITER] = ...
+                           SolveSpectralBIE(axi, bxi, xi, ...
+                                            RotationMatrix, ...
+                                            NGSphere, NGthet, NGphi, ...
+                                            eta, wg, ...
+                                            af, bf, auinf, buinf, ...
+                                            nlat, nlon, mask_a, mask_b, ...
+                                            mu, lam, N, ToleranceGMRES, ...
+                                            cuprev)
+cxi = [axi(mask_a); bxi(mask_b)];
+CHI = reshape(xi,nlat*nlon,3)';
+%% Compute uinf - Gf
+cf = [af(mask_a); bf(mask_b)];
+cGf = computeGf(cf, cxi, CHI, RotationMatrix, ...
+                eta, wg, mask_a, mask_b, mu, ...
+                nlat, nlon, N, NGSphere, NGphi);
+cuinf = [auinf(mask_a); buinf(mask_b)];
+rhs = cuinf - cGf;
+
+%% Compute u + (lam-1)/8/pi*Ku
+Kernel = computeK(cxi, CHI, RotationMatrix, ...
+                  eta, wg, mask_a, mask_b, ...
+                  nlat, nlon, N, NGSphere, NGthet, NGphi);
+
+%% Solve spectal Galerkin BIE 
+[cu, FLAG, RELRES, ITER] = gmres(@computeSpectralBIE,rhs,[], ...
+                                 ToleranceGMRES,3*(N+1)^2,[],[],cuprev);
+
+au = zeros(size(axi));  bu = zeros(size(bxi));
+au(mask_a) = cu(1:3*(N+1)*(N+2)/2);
+bu(mask_b) = cu(3*(N+1)*(N+2)/2+1:3*(N+1)^2);
+u = shsgcm(au, bu);
+    %%
+    function c = computeSpectralBIE(cu)
+        %%
+        au = zeros(size(axi));  bu = zeros(size(bxi));
+        au(mask_a) = cu(1:3*(N+1)*(N+2)/2);
+        bu(mask_b) = cu(3*(N+1)*(N+2)/2+1:3*(N+1)^2);
+        u = shsgcm(au, bu);
+        c = computeKu(Kernel, cu, u, RotationMatrix, ...
+                      mask_a, mask_b, lam, N, nlat, nlon, NGSphere);
+    end
+end