0001 function [u,v,adjointu,adjointv,direction]= ...
0002 grad_J(Rhs,S,F,P,P0,area,specified_velocity,nods,vx_obs,vy_obs,...
0003 index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,B_bar);
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 disp([' Recalculating model velocity with last iteration from velfinder']);
0015 Rhs_parsed=P*(Rhs - F*specified_velocity);
0016 Fprime=P*F*P';
0017
0018 solution=Fprime\Rhs_parsed;
0019 solution=P'*solution + specified_velocity;
0020 u=solution(1:nods);
0021 v=solution(nods+1:2*nods);
0022 disp(' Done');
0023
0024
0025 disp([' Calculating adjoints.']);
0026 Du=[(u-vx_obs).*weighting.*weighting
0027 (v-vy_obs).*weighting.*weighting];
0028
0029 [nu2,nu3]=visc_grad(index,nel,alpha,beta,u,v,B_bar);
0030
0031 ARhs_aprime=-P*(S*Du);
0032
0033 ADprime=P*F*P';
0034
0035 templambda=P'*(ADprime\ARhs_aprime);
0036 adjointu=templambda(1:nods);
0037 adjointv=templambda((nods+1):(nods*2));
0038 disp(' Done');
0039
0040 disp(' Building the gradJ');
0041 value=zeros(nel,27);
0042 row=zeros(nel,27);
0043 count=0;
0044 for m=1:3
0045 for k=1:3
0046
0047 for l=1:3
0048
0049 count=count+1;
0050 row(:,count)=index(:,m);
0051
0052
0053 value(:,count)=-area/3.*z_thick_bar.*nu3.*(...
0054 ( 2.*alpha(:,k).*u(index(:,k)) +beta(:,k).*v(index(:,k)) ) .*2.*alpha(:,l).*adjointu(index(:,l))+...
0055 ( beta(:,k).*u(index(:,k)) + alpha(:,k).*v(index(:,k)) ) .*( beta(:,l).*adjointu(index(:,l)) +alpha(:,l).*adjointv(index(:,l)) )+...
0056 ( 2*beta(:,k).*v(index(:,k)) + alpha(:,k).*u(index(:,k)) ).*2.* beta(:,l).*adjointv(index(:,l))...
0057 );
0058
0059
0060 end
0061 end
0062 end
0063
0064
0065 direction=sparse(row,ones(size(row)),value,nods,1);
0066 direction=full(direction);
0067 disp(' Done');