%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% diffuse2d_opt_imp.m
%
% demonstration: Find optimal importance function J(x) for random walk
%                in 2D domain.  Solving the integral equation by power method.
%
% Wei Cai
% for LLNL CCMS Summer Institute Lectures Aug.1-3, 2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


NGRID = 40; %number of grid points in both directions
r = 0.05;   %absorption rate inside domain
Niter = 1000;

x=[1:NGRID]; y=[1:NGRID];
xx=ones(NGRID,1)*[1:NGRID];
yy=[1:NGRID]'*ones(1,NGRID);

xs=ceil(NGRID*0.25);
ys=ceil(NGRID*0.4);

%absorption function
q=max((abs(xx/NGRID-0.5)>0.42)|(abs(yy/NGRID-0.5)>0.35), r);

%source function
s=zeros(NGRID,NGRID);
s(ys,xs)=1;

%gain function (at domain boundary)
g=(q==1);

%compute new transition probability in importance sampling
dr = [1 0; -1 0; 0 1; 0 -1];
K=ones(NGRID,NGRID,4)/4;   %old transition probability matrix
for i=1:NGRID,
   for j=1:NGRID,
      K(j,i,:)=K(j,i,:)*(1-q(j,i));
   end
end

%solve for psi(x,y)
J=g;
for iter=1:Niter,
   KJ=zeros(size(J));
   KJ1=KJ; KJ1(:,:)=K(:,:,1).*J(:,[2:NGRID,1]);
   KJ2=KJ; KJ2(:,:)=K(:,:,2).*J(:,[NGRID,1:NGRID-1]);
   KJ3=KJ; KJ3(:,:)=K(:,:,3).*J([2:NGRID,1],:);
   KJ4=KJ; KJ4(:,:)=K(:,:,4).*J([NGRID,1:NGRID-1],:);
   KJ=KJ1+KJ2+KJ3+KJ4;
   J=KJ+g;
end

G=sum(sum(s.*J));
disp(sprintf('exact solution G=s*J= %e',G));

%plot results
fs=18;
figure(2);
mesh(log10(J));

set(gca,'FontSize',fs);
xlabel('x');
ylabel('y');
zlabel('J');

