grad_J_flow

PURPOSE ^

GRAD_J - compute the gradient of the misfit with respect to the inversed parameter

SYNOPSIS ^

function [u,v,adjointu,adjointv,direction]=grad_J(Rhs,S,F,P,P0,area,specified_velocity,nods,vx_obs,vy_obs,index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,B_bar);

DESCRIPTION ^

GRAD_J - compute the gradient of the misfit with respect to the inversed parameter

   this routine is calculates the difference between the calculated and the given velocity
   field (cost function), then calculates the adjoint state and uses it to determine the gradient of the
   cost function to update the spatial distribution of viscosity.

   Usage:
      [u,v,adjointu,adjointv,direction]=grad_J(Rhs,S,F,P,P0,area,specified_velocity,nods,vx_obs,vy_obs,index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,B_bar);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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 %GRAD_J - compute the gradient of the misfit with respect to the inversed parameter
0005 %
0006 %   this routine is calculates the difference between the calculated and the given velocity
0007 %   field (cost function), then calculates the adjoint state and uses it to determine the gradient of the
0008 %   cost function to update the spatial distribution of viscosity.
0009 %
0010 %   Usage:
0011 %      [u,v,adjointu,adjointv,direction]=grad_J(Rhs,S,F,P,P0,area,specified_velocity,nods,vx_obs,vy_obs,index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,B_bar);
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 %G_calc;
0031 ARhs_aprime=-P*(S*Du);
0032 %ADprime=P0*(F-G_bar)*P0';
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');

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