objectivefunction_C_flow

PURPOSE ^

OBJECTIVEFUNCTION_C - compute the misfit of a model

SYNOPSIS ^

function [J,u,v]=objectivefunction_C(search_scalar,B,glen_coeff,direction,Rhs,S,F,P,specified_velocity,nods,nel,vx_obs,vy_obs,u_objective,v_objective,index,alpha, beta, area,weighting,rowD,colD,rowDshort,colDshort,valueuu,valuevv,valueuv,matrix_bar,criterion)

DESCRIPTION ^

OBJECTIVEFUNCTION_C - compute the misfit of a model

   this routine updates the value of the spatial distribution of viscosity, calculates the
   new velocity field and compares it to the observed velocity field (J)

   Usage:
      [J,u,v]=objectivefunction_C(search_scalar,B,glen_coeff,direction,Rhs,S,F,P,specified_velocity,nods,nel,vx_obs,vy_obs,u_objective,v_objective,index,alpha, beta, area,weighting,rowD,colD,rowDshort,colDshort,valueuu,valuevv,valueuv,matrix_bar,criterion)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [J,u,v]=objectivefunction_C(search_scalar,B,glen_coeff,direction,Rhs,S,F,P,specified_velocity,...
0002    nods,nel,vx_obs,vy_obs,u_objective,v_objective,index,alpha, beta, area,weighting,rowD,colD,rowDshort,colDshort,valueuu,valuevv,...
0003    valueuv,matrix_bar,criterion)
0004 %OBJECTIVEFUNCTION_C - compute the misfit of a model
0005 %
0006 %   this routine updates the value of the spatial distribution of viscosity, calculates the
0007 %   new velocity field and compares it to the observed velocity field (J)
0008 %
0009 %   Usage:
0010 %      [J,u,v]=objectivefunction_C(search_scalar,B,glen_coeff,direction,Rhs,S,F,P,specified_velocity,nods,nel,vx_obs,vy_obs,u_objective,v_objective,index,alpha, beta, area,weighting,rowD,colD,rowDshort,colDshort,valueuu,valuevv,valueuv,matrix_bar,criterion)
0011 
0012 B_new=B+search_scalar*direction;
0013 %@@@ B_bar=matrix_bar*B_new;
0014 
0015 B_bar=(B_new(index(:,1))+B_new(index(:,2))+B_new(index(:,3)))/3;
0016 %### nu_bar=10^14*ones(nel,1);
0017 
0018 %###AK
0019 nu_bar=viscosity(index,nel,alpha,beta,u_objective,v_objective,B_bar,glen_coeff);
0020 %###
0021 
0022 converged_yet=0;
0023 convergence_count=1;
0024 loop=0;
0025 
0026 
0027 while (~converged_yet)
0028    clear u v
0029    
0030    
0031    nu_bar_uu=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];
0032    nu_bar_vv=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];
0033    nu_bar_uv=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];
0034    
0035    if (loop>100)
0036         u=u_old; v=v_old;
0037       break;
0038    end;
0039    nu_bar_oldvalue=nu_bar;
0040    loop=loop+1;
0041    
0042    Duu=sparse(rowDshort,colDshort,nu_bar_uu.*valueuu,nods,nods);
0043    Duu=Duu+triu(Duu,1)';
0044    
0045    Dvv=sparse(rowDshort,colDshort,nu_bar_vv.*valuevv,nods,nods);
0046    Dvv=Dvv+triu(Dvv,1)';
0047    
0048    F=[Duu sparse(rowD,colD,nu_bar_uv.*valueuv,nods,nods)
0049       sparse(colD,rowD,nu_bar_uv.*valueuv,nods,nods) Dvv]; 
0050    
0051    
0052    Rhs_parsed=P*(Rhs - F*specified_velocity);
0053    Fprime=P*F*P';
0054    
0055    
0056    solution=Fprime\Rhs_parsed;
0057    solution=P'*solution + specified_velocity;
0058    u=solution(1:nods);
0059    v=solution(nods+1:2*nods);
0060    nu_bar=viscosity(index,nel,alpha,beta, u,v,B_bar,glen_coeff);
0061    
0062    %Test for direct shooting convergence
0063    if convergence_count>1,
0064 
0065        dug=[u; v]-[u_old; v_old];
0066        ndug=norm(dug,2);
0067        ug=[u_old;v_old];
0068        nug=norm(ug,2);
0069 
0070        change=ndug/nug;
0071 
0072        %Figure out if viscosity converged
0073        if change>criterion,
0074            converged_yet=0;
0075        else
0076            converged_yet=1;
0077        end
0078    else
0079        converged_yet=0;
0080    end
0081 
0082    u_old=u;
0083    v_old=v;
0084 
0085 convergence_count=convergence_count+1;
0086  
0087 end 
0088 %disp(['# of elements over criterion: ' num2str(sum(test))]);
0089 
0090 %@@@AK save F_file F
0091 save F_file F 
0092 
0093 Du=[(u-vx_obs).*weighting
0094    (v-vy_obs).*weighting];
0095 
0096 
0097 J=Du'*S*Du;

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