1*c4762a1bSJed Brownclear all; 2*c4762a1bSJed Brownclose all; 3*c4762a1bSJed Brown 4*c4762a1bSJed Brownfontsize_labels = 14; 5*c4762a1bSJed Brownfontsize_grid = 12; 6*c4762a1bSJed Brownfontname = 'Times'; 7*c4762a1bSJed Brown 8*c4762a1bSJed Brown%the files optimize** contain 9*c4762a1bSJed Brown%Grad --- gradient, 10*c4762a1bSJed Brown%Init_ts --- initial condition of forward problem 11*c4762a1bSJed Brown%Init_adj--- initial condition of backward problem 12*c4762a1bSJed Brown 13*c4762a1bSJed Brown%additionally 14*c4762a1bSJed Brown%xg - the grid 15*c4762a1bSJed Brown%obj- the objective function 16*c4762a1bSJed Brown%ic - initial condition (the one to be optimized) 17*c4762a1bSJed Brown 18*c4762a1bSJed Brown 19*c4762a1bSJed Brown% run('PDEadjoint/optimize06.m') 20*c4762a1bSJed Brown% figure(2) 21*c4762a1bSJed Brown% plot(xg,obj,'k-','Markersize',10,'Linewidth',2); drawnow 22*c4762a1bSJed Brown% hold on 23*c4762a1bSJed Brown% plot(xg,fwd,'b-','Markersize',10); drawnow 24*c4762a1bSJed Brown% plot(xg,Init_ts,'r*-','Markersize',10); drawnow 25*c4762a1bSJed Brown% plot(xg,Grad,'bs','Markersize',10); drawnow 26*c4762a1bSJed Brown% plot(xg,ic,'g-','LineWidth',2,'Markersize',12); 27*c4762a1bSJed Brown 28*c4762a1bSJed Brown% figure(29) 29*c4762a1bSJed Brown% run('tss.m') 30*c4762a1bSJed Brown% plot(xg,init) 31*c4762a1bSJed Brown% hold on 32*c4762a1bSJed Brown% plot(xg,fin,'ro-') 33*c4762a1bSJed Brown% break 34*c4762a1bSJed Brown% figure(3) 35*c4762a1bSJed Brown% plot(xg,Init_adj,'k*-','Markersize',10); drawnow 36*c4762a1bSJed Brown% hold on 37*c4762a1bSJed Brown% plot(xg,Grad,'go-','Markersize',10); drawnow 38*c4762a1bSJed Brown% plot(xg,Init_adj,'r*-','Markersize',10); drawnow 39*c4762a1bSJed Brown% plot(xg,obj,'r*-','Markersize',10); drawnow 40*c4762a1bSJed Brown% figure(55) 41*c4762a1bSJed Brown% plot(xg,Init_adj,'k*-','Markersize',10); drawnow 42*c4762a1bSJed Brown% hold on 43*c4762a1bSJed Brown% plot(xg,-2*(obj-fwd),'ro','Markersize',10); drawnow 44*c4762a1bSJed Brown 45*c4762a1bSJed Brownfigure(15) 46*c4762a1bSJed Brownfor ii=11:13 47*c4762a1bSJed Brownfile=sprintf('PDEadjoint/optimize%02d.m',ii); 48*c4762a1bSJed Brownrun(file) 49*c4762a1bSJed Brownplot(grid,Init_ts,'ko-','Markersize',10); drawnow; 50*c4762a1bSJed Brownhold on 51*c4762a1bSJed Brown%plot(grid,temp,'c*','Markersize',10); drawnow; 52*c4762a1bSJed Brownplot(grid,exact,'r*','Markersize',10); drawnow; 53*c4762a1bSJed Brown%plot(grid,Curr_sol,'g-','Markersize',10); drawnow; 54*c4762a1bSJed Brown%plot(grid,Init_adj,'bo-','Markersize',10); drawnow; 55*c4762a1bSJed Brownend 56*c4762a1bSJed Brownxlabel('x (GLL grid)'); 57*c4762a1bSJed Brownylabel('f(x)- objective'); 58*c4762a1bSJed Brown% 59*c4762a1bSJed Brown% break 60*c4762a1bSJed Brown% tt=senmask.*obj; 61*c4762a1bSJed Brown% tt(abs(tt)==0)=NaN; 62*c4762a1bSJed Brown% plot(xg,tt,'ks-','Markersize',10); drawnow; 63*c4762a1bSJed Brown% break 64*c4762a1bSJed Brown% Err1=Err; TAO1=TAO; 65*c4762a1bSJed Brown% 66*c4762a1bSJed Brown% %break 67*c4762a1bSJed Brown% figure(15) 68*c4762a1bSJed Brown% for ii=35:40 69*c4762a1bSJed Brown% file=sprintf('PDEadjoint_hc/optimize%02d.m',ii); 70*c4762a1bSJed Brown% run(file) 71*c4762a1bSJed Brown% plot(xg,Init_ts,'go-','Markersize',10); drawnow; 72*c4762a1bSJed Brown% hold on 73*c4762a1bSJed Brown% plot(xg,obj,'r-','Markersize',10); drawnow; 74*c4762a1bSJed Brown% plot(xg,fwd,'bo-','Markersize',10); drawnow; 75*c4762a1bSJed Brown% end 76*c4762a1bSJed Brown% xlabel('x (GLL grid)'); 77*c4762a1bSJed Brown% ylabel('f(x)- objective'); 78*c4762a1bSJed Brown% 79*c4762a1bSJed Brown% 80*c4762a1bSJed Brown% figure(99) 81*c4762a1bSJed Brown% semilogy(Err1,'k-','Markersize',6,'LineWidth',2); drawnow; 82*c4762a1bSJed Brown% hold on 83*c4762a1bSJed Brown% semilogy(Err,'r-','Markersize',6,'LineWidth',2); drawnow; 84*c4762a1bSJed Brown% 85*c4762a1bSJed Brown% 86*c4762a1bSJed Brown% set(gca,'FontName',fontname) 87*c4762a1bSJed Brown% set(gca,'FontSize',fontsize_grid) 88*c4762a1bSJed Brown% set(gca,'FontSize',fontsize_labels) 89*c4762a1bSJed Brown% 90*c4762a1bSJed Brown% legend('Discrete Objective','Continous Objective') 91*c4762a1bSJed Brown% xlabel('Iterations'); 92*c4762a1bSJed Brown% ylabel('Error solution'); 93*c4762a1bSJed Brown% grid on 94*c4762a1bSJed Brown% legend boxoff; 95*c4762a1bSJed Brown% axis tight; axis square 96*c4762a1bSJed Brown% 97*c4762a1bSJed Brown% break 98*c4762a1bSJed Brown% 99*c4762a1bSJed Brown% figure(12) 100*c4762a1bSJed Brown% plot(xg,objk,'b*-','Markersize',6,'LineWidth',2); drawnow; 101*c4762a1bSJed Brown% 102*c4762a1bSJed Brown% hold on 103*c4762a1bSJed Brown% plot(xg,Init_ts,'ro-','Markersize',6,'LineWidth',2); drawnow; 104*c4762a1bSJed Brown% plot(xg,ic,'ks-','Markersize',6,'LineWidth',2); drawnow; 105*c4762a1bSJed Brown% legend('Objective','Optimal','Starting') 106*c4762a1bSJed Brown% legend boxoff 107*c4762a1bSJed Brown% xlabel('GLL grid'); 108*c4762a1bSJed Brown% ylabel('Diffusion solution (Data assimilation)'); 109*c4762a1bSJed Brown% set(gca,'FontName',fontname) 110*c4762a1bSJed Brown% set(gca,'FontSize',fontsize_grid) 111*c4762a1bSJed Brown% set(gca,'FontSize',fontsize_labels) 112*c4762a1bSJed Brown% figure(95) 113*c4762a1bSJed Brown% t=0.6; mu=0.001;x=xg; 114*c4762a1bSJed Brown% plot(xg,2.0*mu*pi*sin(pi*x).*exp(-pi^2*t*mu)./(2.0+exp(-pi^2*t*mu).*cos(pi*x))); 115*c4762a1bSJed Brown% 116*c4762a1bSJed Brown% break 117*c4762a1bSJed Brown% figure(1) 118*c4762a1bSJed Brown% plot(xg,Init_ts,'ro-','Markersize',8,'LineWidth',2); drawnow; 119*c4762a1bSJed Brown% %break 120*c4762a1bSJed Brown% figure(2);set(gca,'FontSize',18);hold on 121*c4762a1bSJed Brown% 122*c4762a1bSJed Brown% run('PDEadjoint/optimize00.m') 123*c4762a1bSJed Brown% plot(xg,Grad,'k*-'); 124*c4762a1bSJed Brown% run('PDEadjoint/optimize04.m') 125*c4762a1bSJed Brown% plot(xg,Grad,'ro-'); 126*c4762a1bSJed Brown% 127*c4762a1bSJed Brown% set(gca,'FontName',fontname) 128*c4762a1bSJed Brown% set(gca,'FontSize',fontsize_grid) 129*c4762a1bSJed Brown% set(gca,'FontSize',fontsize_labels) 130*c4762a1bSJed Brown% 131*c4762a1bSJed Brown% xlabel('x (GLL grid)'); 132*c4762a1bSJed Brown% ylabel('f(x)- objective'); 133*c4762a1bSJed Brown% 134*c4762a1bSJed Brown% legend('Grad at it=0','Grad at it=1') 135*c4762a1bSJed Brown% 136*c4762a1bSJed Brown% figure(10) 137*c4762a1bSJed Brown% run('fd.m') 138*c4762a1bSJed Brown% %plot(gradj) 139*c4762a1bSJed Brown% plot(xg,gradj./Mass,'ro-','LineWidth',2,'Markersize',12); 140*c4762a1bSJed Brown% hold on 141*c4762a1bSJed Brown% run('PDEadjoint/optimize01.m') 142*c4762a1bSJed Brown% plot(xg,Grad,'k*-','LineWidth',2,'Markersize',10); 143*c4762a1bSJed Brown% 144*c4762a1bSJed Brown% set(gca,'FontName',fontname) 145*c4762a1bSJed Brown% set(gca,'FontSize',fontsize_grid) 146*c4762a1bSJed Brown% set(gca,'FontSize',fontsize_labels) 147*c4762a1bSJed Brown% 148*c4762a1bSJed Brown% legend('Gradient FD','Gradient Adjoint') 149*c4762a1bSJed Brown% xlabel('x (GLL grid)'); 150*c4762a1bSJed Brown% ylabel('Gradient'); 151*c4762a1bSJed Brown% axis tight; axis square 152*c4762a1bSJed Brown% 153*c4762a1bSJed Brown% errgrad=max(abs(gradj./Mass-Grad)) 154*c4762a1bSJed Brown% 155*c4762a1bSJed Brown% 156*c4762a1bSJed Brown% figure(21) 157*c4762a1bSJed Brown% semilogy(1:21,TAO,'r','LineWidth',2) 158*c4762a1bSJed Brown% % hold on 159*c4762a1bSJed Brown% % semilogy(1:31,L2,'r','LineWidth',2) 160*c4762a1bSJed Brown% grid on 161*c4762a1bSJed Brown% xlabel('No iterations'); 162*c4762a1bSJed Brown% ylabel('Cost function'); 163*c4762a1bSJed Brown% set(gca,'FontName',fontname) 164*c4762a1bSJed Brown% set(gca,'FontSize',fontsize_grid) 165*c4762a1bSJed Brown% set(gca,'FontSize',fontsize_labels) 166*c4762a1bSJed Brown% % 167*c4762a1bSJed Brown% % legend('TAO','User') 168