cielodiagnostic

PURPOSE ^

diagnostic solution

SYNOPSIS ^

function md=cielodiagnostic(md)

DESCRIPTION ^

diagnostic solution

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function md=cielodiagnostic(md)
0002 %diagnostic solution
0003 
0004     %timing
0005     t1=clock;
0006 
0007     %some configuration transfer between Ice and Cielo
0008     infile='./cielo.dat';
0009 
0010     analysis_dh='diagnostic_horiz';
0011     analysis_dbv='diagnostic_basevert';
0012     analysis_dv='diagnostic_vert';
0013 
0014     %configure, and initialise data model and parameter defaults:
0015     SetParameterDefaults;
0016 
0017     %Build all models requested for diagnostic simulation
0018     %remark, partitions are all identical.
0019     m_dh=CreateFEMModel(md,analysis_dh);
0020     if strcmpi(md.type,'3d'),
0021         m_dbv=CreateFEMModel(md,analysis_dbv);
0022         m_dv=CreateFEMModel(md,analysis_dv);
0023     end
0024 
0025 
0026     %Build partitioning vectors. We only need one partitioning for all our models, we made sure of it in Imp!
0027     part=m_dh.part;
0028     tpart=m_dh.tpart';
0029     indx=1:6:m_dh.uset.gsize; indx=indx(tpart);
0030     indy=2:6:m_dh.uset.gsize; indy=indy(tpart);
0031     indz=3:6:m_dh.uset.gsize; indz=indz(tpart);
0032 
0033     % figure out number of dof: just for information purposes.
0034     md.dof=m_dh.uset.fsize; %biggest dof number
0035 
0036     %rifts
0037     m_dh.rifts=md.rifts;
0038 
0039     %Get horizontal solution.
0040     disp(sprintf('\n%s',['computing horizontal velocities...']));
0041 
0042     u_g=cielodiagnostic_core_nonlinear(m_dh,m_dh.params,{},analysis_dh);
0043 
0044     if strcmpi(md.type,'3d'),
0045     
0046         %swith u_g to workspace partitioning.
0047         u_g=SwitchPartitioning(u_g,'workspace',part,tpart,[1 2]); 
0048         
0049         %extrude velocities for collapsed penta elements
0050         u_g=CieloVelocityExtrude(md,u_g);
0051 
0052         %Compute depth averaged velocity and add it to inputs
0053         velocity_average=CieloHorizontalVelocityDepthAverage(md,u_g);
0054 
0055         %swith u_g and velocity_average to cluster partitioning
0056         u_g=SwitchPartitioning(u_g,'cluster',part,tpart,[1 2]); 
0057         velocity_average=SwitchPartitioning(velocity_average,'cluster',part,tpart,[1 2]); 
0058         
0059 
0060         disp(sprintf('\n%s',['computing basal vertical velocities...']));
0061         SetUset(m_dbv);
0062         u_g_basevert=cielodiagnostic_core_linear(m_dbv,m_dbv.params,struct('velocity',u_g,'velocity_average',velocity_average),analysis_dbv); 
0063         
0064         %use u_g_basevert to constrain vertical velocity
0065         SetUset(m_dv);
0066         m_dv.y_g=u_g_basevert;
0067 
0068         %reduce constraining vector to s-set (the one we solve on)
0069         m_dv.y_s = Reducevectorg( m_dv.y_g);
0070 
0071         %run core linear to solve for vertical velocity
0072         disp(sprintf('\n%s',['computing vertical velocities...']));
0073         u_g_vert=cielodiagnostic_core_linear(m_dv,m_dv.params,struct('velocity',u_g),analysis_dv); 
0074         
0075         %load results onto model: carefull, u_g and u_g_vert are in cluster partitioning
0076         md.vx=u_g(indx)*md.yts;
0077         md.vy=u_g(indy)*md.yts;
0078         md.vz=u_g_vert(indz)*md.yts;
0079         md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
0080         
0081         %Computation of pressure with Pattyn's assumptions (P=rho_ice*g*(s-z) in Pa)
0082         md.pressure=md.rho_ice*md.g*(md.surface-md.z);
0083     else
0084         %load results onto model
0085         md.vx=u_g(indx)*md.yts;
0086         md.vy=u_g(indy)*md.yts;
0087         md.vz=zeros(md.numberofgrids,1);
0088         md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
0089     
0090         %Computation of pressure with Pattyn's assumptions (P=rho_ice*g*(s-z) in Pa)
0091         md.pressure=md.rho_ice*md.g*(md.thickness);
0092     end        
0093 
0094     %timing
0095     t2=clock;
0096     disp(sprintf('\n%s\n',['   solution converged in ' num2str(etime(t2,t1)) ' seconds']));

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