debug

PURPOSE ^

DEBUG - comparison of CIELO and ICE solutions

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

DEBUG - comparison of CIELO and ICE solutions

   This script is used to debug Cielo and Ice by comparing
   the results of both solutions

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 %DEBUG - comparison of CIELO and ICE solutions
0002 %
0003 %   This script is used to debug Cielo and Ice by comparing
0004 %   the results of both solutions
0005 
0006 %first get this into the cielo code, at the core solution level
0007 %PetscSynchronizedPrintf(MPI_COMM_WORLD,"Messing up here\n");
0008 %PetscSynchronizedFlush(MPI_COMM_WORLD);
0009 %BatchDebug(K_gg,p_g,femmodel,"Debug.bin");
0010 %PetscFinalize(); abort();
0011 
0012 
0013 %then watch out for the try, catch statements. Try without them, so see if solve goes through
0014 %nicely, if not correct. When everything works, reinstate the try catch.
0015 
0016 %an alter is used to modify the core solution in Ice, to stop right after Emg. But for cielo,
0017 %we don't have an alter, just use the first 5 lines of this file and put them in the batch
0018 %solution.
0019 
0020 md=model;
0021 md=mesh(md,'DomainOutline.exp',200000);
0022 md=geography(md,'all','');
0023 md=parameterize(md,'Square.par');
0024 md=extrude(md,3,2);
0025 md.eps_rel=.1;
0026 md.eps_abs=NaN;
0027 md=setelementstype(md,'pattyn','Pattyn.exp','fill','macayeal');
0028 %md=setelementstype(md,'fill','macayeal');
0029 md.viscosity_overshoot=0;
0030 md.vx=NaN; md.vy=NaN; md.vz=NaN;
0031 
0032 %first alter solution
0033 !cp /home/larour/Ice/ice1/Solver/Ice/solutions/icediagnostic_core_nonlinear.m ./
0034 !/home/larour/Ice/ice1/Solver/Ice/solutions/DiagnosticAlter.pl
0035 !mv icediagnostic_core_nonlinear.alt.m icediagnostic_core_nonlinear.m
0036 
0037 !rm  -rf systemmatrices_ice.mat
0038 
0039 %Compute solution with Ice model. Catch error triggered to get our of core function
0040 try,
0041     md=solve(md,'diagnostic','ice');
0042 catch
0043 end
0044 
0045 %recover pg and Kgg
0046 results=load('systemmatrices_ice');
0047 Kgg_ice=results.K_gg;
0048 pg_ice=results.p_g;
0049 
0050 %erase alter
0051 %!rm icediagnostic_core_nonlinear.m
0052 ls
0053 
0054 %Now for cielo
0055 md.vx=NaN; md.vy=NaN; md.vz=NaN; md.vel=NaN;
0056 
0057 md.cluster='wilkes';
0058 md.np=3;
0059 md.batch=1;
0060 md.time=10;
0061 md=solve(md,'diagnostic_horiz','cielo');
0062 
0063 cielo_rc_location=which('cielo.rc'); [codepath,executionpath]=ProcessBatchParametersFromCieloRc(md.cluster,cielo_rc_location);
0064 
0065 pause(5);
0066 results=parseresultsfromdisk([executionpath '/Debug.bin']);
0067 
0068 
0069 %recover pg and Kgg
0070 tpart=results{1};
0071 Kgg_cielo=results{2};
0072 pg_cielo=results{3};
0073 
0074 %apply partitioning to Kgg and pg:
0075 gsize=md.numberofgrids*6;
0076 indx=1:6:gsize; indx=indx(tpart);
0077 indy=2:6:gsize; indy=indy(tpart);
0078 indz=3:6:gsize; indz=indz(tpart);
0079 indu=4:6:gsize; indu=indu(tpart);
0080 indv=5:6:gsize; indv=indv(tpart);
0081 indw=6:6:gsize; indw=indw(tpart);
0082 
0083 
0084 ind=zeros(gsize,1); 
0085 ind(1:6:gsize)=indx;
0086 ind(2:6:gsize)=indy;
0087 ind(3:6:gsize)=indz;
0088 ind(4:6:gsize)=indu;
0089 ind(5:6:gsize)=indv;
0090 ind(6:6:gsize)=indw;
0091 
0092 Kgg_cielo=Kgg_cielo(ind,ind);
0093 pg_cielo=pg_cielo(ind);
0094 
0095 disp('Trying to see if stiffness is identical');
0096 
0097 %display some information
0098 disp(sprintf('%s%g%s%g','norm(Kgg_ice): ',norm(Kgg_ice,inf),' norm(Kgg_cielo): ',norm(Kgg_cielo,inf)));
0099 
0100 %find some values that are different
0101 kdiff=(Kgg_ice-Kgg_cielo)./Kgg_ice;
0102 pos=find(isnan(kdiff));kdiff(pos)=0;
0103 pos=find(isinf(kdiff));kdiff(pos)=0;
0104 
0105 pos=find(abs(Kgg_cielo)/norm(Kgg_cielo)>10^-7);
0106 pos2=find(abs(Kgg_cielo(pos))==max(abs(Kgg_cielo(pos))));
0107 pos3=pos(pos2(1));
0108 [i0,j0]=ind2sub(size(Kgg_cielo),pos3);
0109 
0110 disp(['Kgg_ice(' num2str(i0) ',' num2str(j0) ')=' num2str(Kgg_ice(i0,j0)) ' Kgg_cielo(' num2str(i0) ',' num2str(j0) ')=' num2str(Kgg_cielo(i0,j0))]);
0111 
0112 
0113 %figure out grids for these dofs.
0114 grid1=floor((i0-mod(i0,6))/6)+1;
0115 grid2=floor((j0-mod(j0,6))/6)+1;
0116 
0117 %now find elements with grid1 and grid2 in them
0118 [elements1,grids1]=find(md.elements==grid1); 
0119 [pos2,grids2]=find(md.elements(elements1,:)==grid2); 
0120 elements=elements1(pos2);
0121 
0122 disp(['Elements ' num2str(elements') ' Grids ' num2str(grid1) ' ' num2str(grid2) ' Dofs ' num2str(mod(i0,6)) ' ' num2str(mod(j0,6))]);
0123 
0124 
0125 
0126 disp('Trying to see if loads are identical');
0127 disp(sprintf('%s%g%s%g','norm(pg_ice): ',norm(pg_ice),' norm(pg_cielo): ',norm(pg_cielo)));
0128 
0129 pdiff=(pg_ice-pg_cielo)./pg_ice;
0130 pos=find(isnan(pdiff));pdiff(pos)=0;
0131 pos=find(isinf(pdiff));pdiff(pos)=0;
0132 
0133 pos=find(abs(pg_cielo)/norm(pg_cielo)>10^-7);
0134 pos2=find(abs(pdiff(pos))==max(abs(pdiff(pos))));
0135 i0=pos(pos2(1));
0136 disp(['pg_ice(' num2str(i0) ')=' num2str(pg_ice(i0)) ' pg_cielo(' num2str(i0) ')=' num2str(pg_cielo(i0))]);
0137 
0138 %figure out grid for this dofs.
0139 grid1=floor((i0-mod(i0,6))/6)+1;
0140 
0141 %now find elements with grid1
0142 [elements1,grids1]=find(md.elements==grid1); 
0143 
0144 disp(['Elements [' num2str(elements1') '] Grid [' num2str(grid1) '] Dof [' num2str(mod(i0,6)) ']' ]);

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003