0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
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
0029 md.viscosity_overshoot=0;
0030 md.vx=NaN; md.vy=NaN; md.vz=NaN;
0031
0032
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
0040 try,
0041 md=solve(md,'diagnostic','ice');
0042 catch
0043 end
0044
0045
0046 results=load('systemmatrices_ice');
0047 Kgg_ice=results.K_gg;
0048 pg_ice=results.p_g;
0049
0050
0051
0052 ls
0053
0054
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
0070 tpart=results{1};
0071 Kgg_cielo=results{2};
0072 pg_cielo=results{3};
0073
0074
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
0098 disp(sprintf('%s%g%s%g','norm(Kgg_ice): ',norm(Kgg_ice,inf),' norm(Kgg_cielo): ',norm(Kgg_cielo,inf)));
0099
0100
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
0114 grid1=floor((i0-mod(i0,6))/6)+1;
0115 grid2=floor((j0-mod(j0,6))/6)+1;
0116
0117
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
0139 grid1=floor((i0-mod(i0,6))/6)+1;
0140
0141
0142 [elements1,grids1]=find(md.elements==grid1);
0143
0144 disp(['Elements [' num2str(elements1') '] Grid [' num2str(grid1) '] Dof [' num2str(mod(i0,6)) ']' ]);