icetransient2d

PURPOSE ^

ICETRANSIENT2D - transient solution of 2d models

SYNOPSIS ^

function md=icetransient2d(md);

DESCRIPTION ^

ICETRANSIENT2D - transient solution of 2d models

   Usage:
      md=icetransient2d(md)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function md=icetransient2d(md);
0002 %ICETRANSIENT2D - transient solution of 2d models
0003 %
0004 %   Usage:
0005 %      md=icetransient2d(md)
0006 
0007 global gridset
0008 
0009 %Create fem structure (input of icediagnostic2d)
0010 fem=struct();
0011 
0012 %Figure out which type of elements are present
0013 [fem.ishutter,fem.ismacayealpattyn,fem.isstokes]=DiagnosticSolutionType(md.elements_type);
0014 
0015 %First, build elements,grids,loads, etc ... for horizontal model
0016 if fem.ishutter,
0017     fem.m_ss=CreateFemModel(md,'surface_slope_compute');
0018     fem.m_dhu=CreateFemModel(md,'diagnostic_hutter');
0019 end
0020 if fem.ismacayealpattyn,
0021     fem.m_dh=CreateFemModel(md,'diagnostic_horiz');
0022 end
0023 m_p=CreateFemModel(md,'prognostic');
0024 
0025 %initialize (velocity,pressure)
0026 solution.u_g=zeros(m_p.gridset.gsize,1);
0027 if ~isnan(md.vx), solution.u_g(1:6:m_p.gridset.gsize)=md.vx/md.yts; end
0028 if ~isnan(md.vy), solution.u_g(2:6:m_p.gridset.gsize)=md.vy/md.yts; end
0029 
0030 %initialize geometry and time
0031 solution.thickness=zeros(m_p.gridset.gsize,1);solution.thickness(1:6:m_p.gridset.gsize)=md.thickness;
0032 solution.surface=zeros(m_p.gridset.gsize,1);solution.surface(1:6:m_p.gridset.gsize)=md.surface;
0033 solution.bed=zeros(m_p.gridset.gsize,1);solution.bed(1:6:m_p.gridset.gsize)=md.bed;
0034 solution.time=0;
0035 
0036 %melting and accumulation are taken constant for all iterations
0037 melting=zeros(m_p.gridset.gsize,1);
0038 accumulation=zeros(m_p.gridset.gsize,1);
0039 melting(1:6:end)=md.melting/md.yts;           %from m/yr to m/s
0040 accumulation(1:6:end)=md.accumulation/md.yts; %from m/yr to m/s
0041 
0042 %build constant inputs
0043 fem.inputs=struct();
0044 fem.inputs.accumulation=accumulation;
0045 fem.inputs.melting=melting;
0046 fem.inputs.dt=md.dt;
0047 
0048 %first time step is given by model.
0049 dt=md.dt;
0050 finaltime=md.ndt;
0051 time=dt;
0052 n=1; %counter
0053 
0054 while  time<finaltime+dt, %make sure we run up to finaltime.
0055 
0056     disp(sprintf('\n%s%g%s%g%s%g\n','time [yr]: ',time/md.yts,'    iteration number: ',n,'/',floor(finaltime/dt)));
0057     
0058     solution(n+1).time=time;
0059 
0060     %update inputs
0061     fem.inputs.velocity=solution(n).u_g;
0062     fem.inputs.thickness=solution(n).thickness;
0063     fem.inputs.bed=solution(n).bed;
0064     fem.inputs.surface=solution(n).surface;
0065 
0066     %Deal with velocities.
0067 
0068     %Get horizontal solution.
0069     solution(n+1).u_g=icediagnostic2d(md,fem);
0070 
0071     %compute new thickness
0072     disp(sprintf('%s','   computing new thickness...'));
0073     fem.inputs.velocity_average=solution(n+1).u_g;
0074     new_thickness=iceprognostic_core(m_p,'prognostic',fem.inputs);
0075 
0076     %update surface and bed using the new thickness
0077     disp(sprintf('%s','   updating geometry...'));
0078     indx=1:6:m_p.gridset.gsize; indx=indx(m_p.tpart);
0079     [new_bed,new_surface,new_thickness]=UpdateGeometry(md,new_thickness(indx),solution(n).bed(indx),solution(n).surface(indx));
0080     
0081     %Record bed surface and thickness in the solution
0082     solution(n+1).thickness=zeros(m_p.gridset.gsize,1);solution(n+1).thickness(1:6:m_p.gridset.gsize)=new_thickness;
0083     solution(n+1).bed=zeros(m_p.gridset.gsize,1);solution(n+1).bed(1:6:m_p.gridset.gsize)=new_bed;
0084     solution(n+1).surface=zeros(m_p.gridset.gsize,1);solution(n+1).surface(1:6:m_p.gridset.gsize)=new_surface;
0085 
0086     %figure out if time stepping is good
0087     disp(sprintf('%s','   checking time stepping...'));
0088     [back,dt,time]=TimeStepping(md,solution,dt,time);
0089     if back,
0090         continue;
0091     end
0092 
0093     %update time and counter
0094     time=time+dt;
0095     n=n+1;
0096 
0097 end
0098 
0099 %load results onto model
0100 indx=1:6:m_p.gridset.gsize; indx=indx(m_p.tpart);
0101 indy=2:6:m_p.gridset.gsize; indy=indy(m_p.tpart);
0102 for i=1:length(solution),
0103     solution2(i).vx=solution(i).u_g(indx)*md.yts;
0104     solution2(i).vy=solution(i).u_g(indy)*md.yts;
0105     solution2(i).vel=sqrt(solution2(i).vx.^2+solution2(i).vy.^2);
0106     solution2(i).thickness=solution(i).thickness(indx);
0107     solution2(i).surface=solution(i).surface(indx);
0108     solution2(i).bed=solution(i).bed(indx);
0109     solution2(i).time=solution(i).time/md.yts;
0110 end
0111 md.transient_results=solution2;

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