Index: sm/trunk-jpl/src/m/contrib/larour/tres.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/larour/tres.m (revision 21061) +++ (revision ) @@ -1,115 +1,0 @@ -function md=tres(md,string) -%TRES - transfer results results to corresponding model fields. -% -% Usage: md=tres(md,string) -% -% Example: md=tres(md,'stressbalance'); - -%check number of arguments - -if strcmpi(string,'stressbalance'), - if strcmp(domaintype(md.mesh),'2Dhorizontal'), - md.initialization.vx=md.results.StressbalanceSolution.Vx; - md.initialization.vy=md.results.StressbalanceSolution.Vy; - else - md.initialization.vx=md.results.StressbalanceSolution.Vx; - md.initialization.vy=md.results.StressbalanceSolution.Vy; - md.initialization.vz=md.results.StressbalanceSolution.Vz; - end - md.initialization.vel=md.results.StressbalanceSolution.Vel; - - if isfield(md.results.StressbalanceSolution,'Pressure'), - md.initialization.pressure=md.results.StressbalanceSolution.Pressure; - end - if ~isempty(md.rifts.riftstruct), - if isfield(md.results.StressbalanceSolution,'riftproperties'), - md.rifts.riftproperties=md.results.StressbalanceSolution.riftproperties; - end - end - -elseif strcmpi(string,'dakota'), - md.qmu.results=md.results.dakota; - -elseif strcmpi(string,'flaim'), - md.flaim.solution=md.results.FlaimSolution.solution; - md.flaim.quality =md.results.FlaimSolution.quality; - -elseif strcmpi(string,'transient'), - results=md.results.TransientSolution; - results2.Vel=NaN; - count=1; - for i=1:length(results), - if ~isempty(md.results.TransientSolution(i).Vel), - results2(count).Vel=md.results.TransientSolution(i).Vel; - results2(count).Surface=md.results.TransientSolution(i).Surface; - results2(count).Thickness=md.results.TransientSolution(i).Thickness; - results2(count).Bed=md.results.TransientSolution(i).Bed; - results2(count).Vx=md.results.TransientSolution(i).Vx; - results2(count).Vy=md.results.TransientSolution(i).Vy; - results2(count).time=md.results.TransientSolution(i).time; - results2(count).step=md.results.TransientSolution(i).step; - if ~strcmpi(md.groundingline.migration,'None'), - results2(count).ElementOnIceShelf=md.results.TransientSolution(i).ElementOnIceShelf; - end - count=count+1; - end - end - md.results.TransientSolution=results2; - clear results,results2; -elseif strcmpi(string,'steadystate'), - md.initialization.vx=md.results.SteadystateSolution.Vx; - md.initialization.vy=md.results.SteadystateSolution.Vy; - if isfield(md.results.SteadystateSolution,'Vz'), - md.initialization.vz=md.results.SteadystateSolution.Vz; - end - - md.initialization.vel=md.results.SteadystateSolution.Vel; - md.initialization.pressure=md.results.SteadystateSolution.Pressure; - md.initialization.temperature=md.results.SteadystateSolution.Temperature; - md.basalforcings.groundedice_melting_rate=md.results.SteadystateSolution.BasalforcingsGroundediceMeltingRate; - - if md.inversion.iscontrol==1, - for control_parameters=md.inversion.control_parameters - md.(EnumToModelField(control_parameters))=md.results.SteadystateSolution.(EnumToString(control_parameters)); - end - end - -elseif strcmpi(string,'thermal'), - md.initialization.temperature=md.results.ThermalSolution.Temperature; - md.basalforcings.groundedice_melting_rate=md.results.ThermalSolution.BasalGroundediceMeltingRate; -elseif strcmpi(string,'hydrology'), - md.initialization.watercolumn=md.results.HydrologySolution.Watercolumn; - -else - error(['tres error message: analysis ' string ' not supported yet!']); -end -end -function string=EnumToModelField(enum) % {{{ - %ENUMTOMODELFIELD - output string of model field associated to enum - % - % Usage: - % string=EnumToModelField(enum) - - disp('Warning: EnumToModelField is deprecated, it cannot work with new model definition. This function will be removed in the future'); - - switch enum, - - case ThicknessEnum(), string='thickness'; return - case FrictionCoefficientEnum(), string='drag_coefficient'; return - case MaterialsRheologyBEnum(), string='rheology_B'; return - case MaterialsRheologyBbarEnum(), string='rheology_B'; return - case MaterialsRheologyZEnum(), string='rheology_Z'; return - case MaterialsRheologyZbarEnum(), string='rheology_Z'; return - case BalancethicknessThickeningRateEnum(), string='dhdt'; return - case VxEnum(), string='vx'; return - case InversionVxObsEnum(), string='vx_obs'; return - case VyEnum(), string='vy'; return - case InversionVyObsEnum(), string='vy_obs'; return - case BasalforcingsMeltingRateEnum(), string='basal_melting_rate'; return - case SmbAccumulationRateEnum(), string='surface_accumulation_rate'; return - case SmbAblationRateEnum(), string='surface_ablation_rate'; return - case SmbMassBalanceEnum(), string='surface_mass_balance'; return - otherwise, error(['Enum ' num2str(enum) ' not found associated to any model field']); - - end - end % }}} Index: /issm/trunk-jpl/src/m/solve/outbinread.m =================================================================== --- /issm/trunk-jpl/src/m/solve/outbinread.m (revision 21062) +++ /issm/trunk-jpl/src/m/solve/outbinread.m (revision 21062) @@ -0,0 +1,161 @@ +function results = outbinread(md,filename,steps), +%OUTBINREAD - read specific steps in outbin file +% +% Usage: +% results = outbinread(md,filename,steps) +% +% Example: +% results = outbinread(md,'dir/test.outbin',1:10) + + +if ~exist(filename,'file'), + error(['Binary file ' filename ' not found']); +end + +%Initialize output +results=struct(); + +%Open file +fid=fopen(filename,'rb'); +if(fid==-1), + error(['could not open ',filename,' for binary reading']); +end + +%Read fields until the end of the file. +result = ReadData(fid,md); +if isempty(result), error(['no results found in binary file ' filename]); end +check_nomoresteps=0; +counter = 1; +step = result.step; +while ~isempty(result), + + if check_nomoresteps, + %check that the new result does not add a step, which would be an error: + if result.step>=1, + error('parsing results for a steady-state core, which incorporates transient results!'); + end + end + + %Check step, increase counter if this is a new step + if(step~=result.step & result.step>1) + counter = counter + 1; + step = result.step; + end + + %Add result + if(result.step==0), + %if we have a step = 0, this is a steady state solution, don't expect more steps. + index = 1; + check_nomoresteps=1; + elseif(result.step==1), + index = 1; + else + index = counter; + end + pos=find(result.step==steps); + if ~isempty(pos), + results(pos).(result.fieldname)=result.field; + if(result.step~=-9999), + results(pos).step=result.step; + end + if(result.time~=-9999), + results(pos).time=result.time; + end + end + if result.step>max(steps(:)), + %we are done! + return; + end + + %read next result + try, + result = ReadData(fid,md); + catch me, + disp('WARNING: file corrupted, trying partial recovery'); + result=[]; + end + +end + +fclose(fid); + +function result=ReadData(fid,md) % {{{ + +%read field +[length,count]=fread(fid,1,'int'); + +if count==0, + result=struct([]); +else + fieldname=fread(fid,length,'char'); + fieldname=fieldname(1:end-1)'; + fieldname=char(fieldname); + time=fread(fid,1,'double'); + step=fread(fid,1,'int'); + + type=fread(fid,1,'int'); + M=fread(fid,1,'int'); + if type==1, + field=fread(fid,M,'double'); + elseif type==2, + field=fread(fid,M,'char'); + field=char(field(1:end-1)'); + elseif type==3, + N=fread(fid,1,'int'); + field=transpose(fread(fid,[N M],'double')); + else + error(['cannot read data of type ' num2str(type) ]); + end + + %Process units here FIXME: this should not be done here! + yts=md.constants.yts; + if strcmp(fieldname,'BalancethicknessThickeningRate'), + field = field*yts; + elseif strcmp(fieldname,'HydrologyWaterVx'), + field = field*yts; + elseif strcmp(fieldname,'HydrologyWaterVy'), + field = field*yts; + elseif strcmp(fieldname,'Vx'), + field = field*yts; + elseif strcmp(fieldname,'Vy'), + field = field*yts; + elseif strcmp(fieldname,'Vz'), + field = field*yts; + elseif strcmp(fieldname,'Vel'), + field = field*yts; + elseif strcmp(fieldname,'BasalforcingsGroundediceMeltingRate'), + field = field*yts; + elseif strcmp(fieldname,'BasalforcingsFloatingiceMeltingRate'), + field = field*yts; + elseif strcmp(fieldname,'TotalFloatingBmb'), + field = field/10.^12*yts; %(GigaTon/year) + elseif strcmp(fieldname,'TotalGroundedBmb'), + field = field/10.^12*yts; %(GigaTon/year) + elseif strcmp(fieldname,'TotalSmb'), + field = field/10.^12*yts; %(GigaTon/year) + elseif strcmp(fieldname,'SmbMassBalance'), + field = field*yts; + elseif strcmp(fieldname,'SmbPrecipitation'), + field = field*yts; + elseif strcmp(fieldname,'SmbRunoff'), + field = field*yts; + elseif strcmp(fieldname,'SmbCondensation'), + field = field*yts; + elseif strcmp(fieldname,'SmbAccumulation'), + field = field*yts; + elseif strcmp(fieldname,'SmbMelt'), + field = field*yts; + elseif strcmp(fieldname,'CalvingCalvingrate'), + field = field*yts; + end + + result.fieldname=fieldname; + result.time=time; + if result.time~=-9999, + result.time=time/yts; + end + result.step=step; + result.field=field; + +end +% }}}