iceparameters

PURPOSE ^

ICEPARAMETERS - compute parameters

SYNOPSIS ^

function md=iceparameters(md);

DESCRIPTION ^

ICEPARAMETERS - compute parameters

   This routine is caled by the solver to compute the parameters
   that are given in md.outputparameters (stress, strain rate,...)

   Usage:
      md=iceparameters(md)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function md=iceparameters(md);
0002 %ICEPARAMETERS - compute parameters
0003 %
0004 %   This routine is caled by the solver to compute the parameters
0005 %   that are given in md.outputparameters (stress, strain rate,...)
0006 %
0007 %   Usage:
0008 %      md=iceparameters(md)
0009 
0010 %define global variables
0011 iceglobal
0012 
0013 %determine if run is parallel
0014 if strcmpi(md.cluster,'yes'), cluster=1; else cluster=0;end;
0015 
0016 %for now, only serial support is in
0017 if cluster,
0018     error('iceparameters error message: parallel support not implemented yet');
0019 end
0020 
0021 if strcmpi(md.type,'2d'),
0022     %First, build elements,grids,loads, etc ... for thermal and melting
0023     m_dh=CreateFemModel(md,'diagnostic_horiz');
0024 
0025     %Are any more computations requested?
0026     if ~isempty(md.parameteroutput),
0027         disp(sprintf('\n%s',['computing outputs...']));
0028         %velocity field
0029         u_g=zeros(gridset.gsize,1); 
0030         u_g(1:6:gridset.gsize)=md.vx/md.yts; %from m/a to m/s
0031         u_g(2:6:gridset.gsize)=md.vy/md.yts;
0032 
0033         %Go through fields of parameteroutput and compute
0034         for i=1:length(md.parameteroutput),
0035             if strcmpi(md.parameteroutput{i},'viscousheating'),
0036                 disp(sprintf('%s',['   viscousheating']));
0037                 md.viscousheating=ViscousHeatingCompute(m_dh,struct('velocity',u_g),md.type);
0038             elseif strcmpi(md.parameteroutput{i},'pressure_elem'),
0039                 disp(sprintf('%s',['   pressure_elem']));
0040                 md.pressure_elem=PressureElemCompute(m_dh,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
0041             elseif strcmpi(md.parameteroutput{i},'strainrate'),
0042                 disp(sprintf('%s',['   strainrate']));
0043                 md.strainrate=StrainRateCompute(m_dh,struct('velocity',u_g),md.type);
0044             elseif strcmpi(md.parameteroutput{i},'stress'),
0045                 disp(sprintf('%s',['   stress']));
0046                 md.stress=StressCompute(m_dh,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
0047             elseif strcmpi(md.parameteroutput{i},'stress_bed'),
0048                 disp(sprintf('%s',['   stress_bed']));
0049                 md.stress_bed=StressBedCompute(m_dh,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
0050             elseif strcmpi(md.parameteroutput{i},'stress_surface'),
0051                 disp(sprintf('%s',['   stress_surface']));
0052                 md.stress_surface=StressSurfaceCompute(m_dh,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
0053             elseif strcmpi(md.parameteroutput{i},'deviatoricstress'),
0054                 disp(sprintf('%s',['   deviatoricstress']));
0055                 md.deviatoricstress=DeviatoricStressCompute(m_dh,struct('velocity',u_g),md.type);
0056             end
0057         end
0058     end
0059 
0060 else
0061 
0062     %First, build elements,grids,loads, etc ... for horizontal, base vertical and vertical model
0063     m_dv=CreateFemModel(md,'diagnostic_vert');
0064 
0065     %Are any more computations requested?
0066     if ~isempty(md.parameteroutput),
0067         disp(sprintf('\n%s',['computing outputs...']));
0068         %velocity field
0069         u_g=zeros(gridset.gsize,1); 
0070         u_g(1:6:gridset.gsize)=md.vx/md.yts; %from m/a to m/s
0071         u_g(2:6:gridset.gsize)=md.vy/md.yts;
0072         u_g(3:6:gridset.gsize)=md.vz/md.yts;
0073         if ~isnan(md.pressure),
0074             u_g(4:6:gridset.gsize)=md.pressure;
0075         end
0076 
0077         %Go through fields of parameteroutput and compute
0078         for i=1:length(md.parameteroutput),
0079             if strcmpi(md.parameteroutput{i},'viscousheating'),
0080                 disp(sprintf('%s',['   viscousheating']));
0081                 md.viscousheating=ViscousHeatingCompute(m_dv,struct('velocity',u_g),md.type);
0082             elseif strcmpi(md.parameteroutput{i},'pressure_elem'),
0083                 disp(sprintf('%s',['   pressure_elem']));
0084                 md.pressure_elem=PressureElemCompute(m_dv,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
0085             elseif strcmpi(md.parameteroutput{i},'strainrate'),
0086                 disp(sprintf('%s',['   strainrate']));
0087                 md.strainrate=StrainRateCompute(m_dv,struct('velocity',u_g),md.type);
0088             elseif strcmpi(md.parameteroutput{i},'stress'),
0089                 disp(sprintf('%s',['   stress']));
0090                 md.stress=StressCompute(m_dv,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
0091             elseif strcmpi(md.parameteroutput{i},'stress_bed'),
0092                 disp(sprintf('%s',['   stress_bed']));
0093                 md.stress_bed=StressBedCompute(m_dv,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
0094             elseif strcmpi(md.parameteroutput{i},'stress_surface'),
0095                 disp(sprintf('%s',['   stress_surface']));
0096                 md.stress_surface=StressSurfaceCompute(m_dv,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
0097             elseif strcmpi(md.parameteroutput{i},'deviatoricstress'),
0098                 disp(sprintf('%s',['   deviatoricstress']));
0099                 md.deviatoricstress=DeviatoricStressCompute(m_dv,struct('velocity',u_g),md.type);
0100             end
0101         end
0102     end
0103 
0104 end

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