%DEBUG - comparison of CIELO and ICE solutions
%
%   This script is used to debug Cielo and Ice by comparing
%   the results of both solutions

%first get this into the cielo code, at the core solution level
%PetscSynchronizedPrintf(MPI_COMM_WORLD,"Messing up here\n");
%PetscSynchronizedFlush(MPI_COMM_WORLD);
%BatchDebug(K_gg,p_g,femmodel,"Debug.bin");
%PetscFinalize(); abort();


%then watch out for the try, catch statements. Try without them, so see if solve goes through 
%nicely, if not correct. When everything works, reinstate the try catch.

%an alter is used to modify the core solution in Ice, to stop right after Emg. But for cielo, 
%we don't have an alter, just use the first 5 lines of this file and put them in the batch 
%solution.

md=model;
md=mesh(md,'DomainOutline.exp',200000);
md=geography(md,'all','');
md=parameterize(md,'Square.par');
md=extrude(md,3,2);
md.eps_rel=.1;
md.eps_abs=NaN;
md=setelementstype(md,'pattyn','Pattyn.exp','fill','macayeal');
%md=setelementstype(md,'fill','macayeal');
md.viscosity_overshoot=0;
md.vx=NaN; md.vy=NaN; md.vz=NaN;

%first alter solution 
!cp /home/larour/Ice/ice1/Solver/Ice/solutions/icediagnostic_core_nonlinear.m ./
!/home/larour/Ice/ice1/Solver/Ice/solutions/DiagnosticAlter.pl
!mv icediagnostic_core_nonlinear.alt.m icediagnostic_core_nonlinear.m

!rm  -rf systemmatrices_ice.mat

%Compute solution with Ice model. Catch error triggered to get our of core function
try,
	md=solve(md,'diagnostic','ice');
catch
end

%recover pg and Kgg
results=load('systemmatrices_ice');
Kgg_ice=results.K_gg;
pg_ice=results.p_g;

%erase alter
%!rm icediagnostic_core_nonlinear.m 
ls

%Now for cielo
md.vx=NaN; md.vy=NaN; md.vz=NaN; md.vel=NaN;

md.cluster='wilkes';
md.np=3;
md.batch=1;
md.time=10;
md=solve(md,'diagnostic_horiz','cielo');

cielo_rc_location=which('cielo.rc'); [codepath,executionpath]=ProcessBatchParametersFromCieloRc(md.cluster,cielo_rc_location);

pause(5);
results=parseresultsfromdisk([executionpath '/Debug.bin']);


%recover pg and Kgg
tpart=results{1};
Kgg_cielo=results{2};
pg_cielo=results{3};

%apply partitioning to Kgg and pg:
gsize=md.numberofgrids*6;
indx=1:6:gsize; indx=indx(tpart);
indy=2:6:gsize; indy=indy(tpart);
indz=3:6:gsize; indz=indz(tpart);
indu=4:6:gsize; indu=indu(tpart);
indv=5:6:gsize; indv=indv(tpart);
indw=6:6:gsize; indw=indw(tpart);


ind=zeros(gsize,1); 
ind(1:6:gsize)=indx;
ind(2:6:gsize)=indy;
ind(3:6:gsize)=indz;
ind(4:6:gsize)=indu;
ind(5:6:gsize)=indv;
ind(6:6:gsize)=indw;

Kgg_cielo=Kgg_cielo(ind,ind);
pg_cielo=pg_cielo(ind);

disp('Trying to see if stiffness is identical');

%display some information
disp(sprintf('%s%g%s%g','norm(Kgg_ice): ',norm(Kgg_ice,inf),' norm(Kgg_cielo): ',norm(Kgg_cielo,inf)));

%find some values that are different
kdiff=(Kgg_ice-Kgg_cielo)./Kgg_ice;
pos=find(isnan(kdiff));kdiff(pos)=0;
pos=find(isinf(kdiff));kdiff(pos)=0;

pos=find(abs(Kgg_cielo)/norm(Kgg_cielo)>10^-7);
pos2=find(abs(Kgg_cielo(pos))==max(abs(Kgg_cielo(pos))));
pos3=pos(pos2(1));
[i0,j0]=ind2sub(size(Kgg_cielo),pos3);

disp(['Kgg_ice(' num2str(i0) ',' num2str(j0) ')=' num2str(Kgg_ice(i0,j0)) ' Kgg_cielo(' num2str(i0) ',' num2str(j0) ')=' num2str(Kgg_cielo(i0,j0))]);


%figure out grids for these dofs.
grid1=floor((i0-mod(i0,6))/6)+1;
grid2=floor((j0-mod(j0,6))/6)+1;

%now find elements with grid1 and grid2 in them
[elements1,grids1]=find(md.elements==grid1); 
[pos2,grids2]=find(md.elements(elements1,:)==grid2); 
elements=elements1(pos2);

disp(['Elements ' num2str(elements') ' Grids ' num2str(grid1) ' ' num2str(grid2) ' Dofs ' num2str(mod(i0,6)) ' ' num2str(mod(j0,6))]);



disp('Trying to see if loads are identical');
disp(sprintf('%s%g%s%g','norm(pg_ice): ',norm(pg_ice),' norm(pg_cielo): ',norm(pg_cielo)));

pdiff=(pg_ice-pg_cielo)./pg_ice;
pos=find(isnan(pdiff));pdiff(pos)=0;
pos=find(isinf(pdiff));pdiff(pos)=0;

pos=find(abs(pg_cielo)/norm(pg_cielo)>10^-7);
pos2=find(abs(pdiff(pos))==max(abs(pdiff(pos))));
i0=pos(pos2(1));
disp(['pg_ice(' num2str(i0) ')=' num2str(pg_ice(i0)) ' pg_cielo(' num2str(i0) ')=' num2str(pg_cielo(i0))]);

%figure out grid for this dofs.
grid1=floor((i0-mod(i0,6))/6)+1;

%now find elements with grid1 
[elements1,grids1]=find(md.elements==grid1); 

disp(['Elements [' num2str(elements1') '] Grid [' num2str(grid1) '] Dof [' num2str(mod(i0,6)) ']' ]);
