icetransient3d

PURPOSE ^

ICETRANSIENT3D - transient solution of 3d models

SYNOPSIS ^

function md=icetransient3d(md);

DESCRIPTION ^

ICETRANSIENT3D - transient solution of 3d models

   Usage:
      md=icetransient3d(md)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function md=icetransient3d(md);
0002 %ICETRANSIENT3D - transient solution of 3d models
0003 %
0004 %   Usage:
0005 %      md=icetransient3d(md)
0006 
0007 global gridset
0008 
0009 %Create fem structure
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, base vertical and vertical model
0016 if fem.ismacayealpattyn,
0017     fem.m_dh=CreateFemModel(md,'diagnostic_horiz');
0018 end
0019 if fem.ishutter,
0020     fem.m_ss=CreateFemModel(md,'surface_slope_compute');
0021     fem.m_dhu=CreateFemModel(md,'diagnostic_hutter');
0022 end
0023 if fem.isstokes,
0024     fem.m_bs=CreateFemModel(md,'bed_slope_compute');
0025     fem.m_ds=CreateFemModel(md,'diagnostic_stokes');
0026 end
0027 fem.m_dbv=CreateFemModel(md,'diagnostic_basevert');
0028 fem.m_dv=CreateFemModel(md,'diagnostic_vert');
0029 
0030 %Second, build elements,grids,loads, etc ... for thermal and prognostic model
0031 fem.m_t=CreateFemModel(md,'thermaltransient');
0032 fem.m_m=CreateFemModel(md,'meltingtransient');
0033 fem.m_p=CreateFemModel(md,'prognostic');
0034 
0035 %initialize (velocity,pressure,...)
0036 solution.u_g=zeros(fem.m_dhu.gridset.gsize,1);
0037 if ~isnan(md.vx), solution.u_g(1:6:fem.m_dhu.gridset.gsize)=md.vx/md.yts; end
0038 if ~isnan(md.vy), solution.u_g(2:6:fem.m_dhu.gridset.gsize)=md.vy/md.yts; end
0039 if ~isnan(md.vz), solution.u_g(3:6:fem.m_dhu.gridset.gsize)=md.vz/md.yts; end
0040 solution.pressure=zeros(fem.m_dhu.gridset.gsize,1); solution.pressure(1:6:fem.m_dhu.gridset.gsize)=md.rho_ice*md.g*(md.surface-md.z); %lithostatic pressure to spin up
0041 solution.t_g=zeros(fem.m_t.gridset.gsize,1);
0042 if ~isempty(md.temperature), 
0043     solution.t_g(1:6:fem.m_t.gridset.gsize)=md.temperature;
0044 else
0045     solution.t_g(1:6:fem.m_t.gridset.gsize)=md.meltingpoint;
0046 end
0047 pos=find(isnan(solution.t_g));
0048 solution.t_g(pos)=md.meltingpoint;
0049 if ~isempty(md.melting), solution.melting_g(1:6:fem.m_m.gridset.gsize)=md.melting;end;
0050 
0051 %initialize geometry and time
0052 solution.mesh=struct('elements',md.elements,'grids',[md.x md.y md.z]);
0053 solution.thickness=zeros(fem.m_dhu.gridset.gsize,1);solution.thickness(1:6:fem.m_dhu.gridset.gsize)=md.thickness;
0054 solution.surface=zeros(fem.m_dhu.gridset.gsize,1);solution.surface(1:6:fem.m_dhu.gridset.gsize)=md.surface;
0055 solution.bed=zeros(fem.m_dhu.gridset.gsize,1);solution.bed(1:6:fem.m_dhu.gridset.gsize)=md.bed;
0056 solution.time=0;
0057 
0058 %build constant inputs
0059 fem.inputs=struct();
0060 
0061 fem.inputs.melting=zeros(fem.m_m.gridset.gsize,1); melting=ShiftLayers(md,md.melting,1,1); fem.inputs.melting(1:6:fem.m_m.gridset.gsize)=melting;
0062 fem.inputs.accumulation=zeros(fem.m_m.gridset.gsize,1); accumulation=ShiftLayers(md,md.accumulation,md.numlayers,1); fem.inputs.accumulation(1:6:fem.m_m.gridset.gsize)=accumulation;
0063 fem.inputs.dt=md.dt;
0064 
0065 %first time step is given by model.
0066 dt=md.dt;
0067 finaltime=md.ndt;
0068 time=dt;
0069 n=1; %counter
0070 
0071 while  time<finaltime+dt, %make sure we run up to finaltime.
0072 
0073     disp(sprintf('\n%s%g%s%g%s%g\n','time [yr]: ',time/md.yts,'    iteration number: ',n,'/',floor(finaltime/dt)));
0074     
0075     solution(n+1).time=time;
0076 
0077     %update inputs
0078     fem.inputs.velocity=solution(n).u_g;
0079     fem.inputs.pressure=solution(n).pressure;
0080     fem.inputs.temperature=full(solution(n).t_g);
0081     fem.inputs.thickness=solution(n).thickness;
0082     fem.inputs.bed=solution(n).bed;
0083     fem.inputs.surface=solution(n).surface;
0084 
0085     %%Deal with temperature first
0086     %disp('   computing temperature...');
0087     %[solution(n+1).t_g loads_t melting_offset]=icethermal_core(fem.m_t,'thermaltransient',fem.inputs);
0088     solution(n+1).t_g=solution(1).t_g;
0089     %disp('   computing melting...');
0090     %fem.inputs.temperature=full(solution(n+1).t_g);
0091     %fem.inputs.melting_offset=melting_offset;
0092     %solution(n+1).melting_g=icemelting_core(fem.m_m,'thermaltransient',fem.inputs);
0093     solution(n+1).melting_g=solution(1).melting_g;
0094 
0095     %Compute depth averaged temperature for 2d type elements.
0096     %temperature_average2d=DepthAverage(md,solution(n+1).t_g(1:6:fem.m_t.gridset.gsize));
0097     %temperature_average=project3d(md,temperature_average2d,'node');
0098     %temperature_average_g=zeros(fem.m_t.gridset.gsize,1);temperature_average_g(1:6:fem.m_t.gridset.gsize)=temperature_average;
0099     %fem.inputs.temperature_average=temperature_average_g;
0100     
0101     %Deal with velocities.
0102     disp(sprintf('%s','   computing velocities...'));
0103     solution(n+1).u_g=icediagnostic3d(md,fem);
0104     solution(n+1).pressure=zeros(fem.m_dv.gridset.gsize,1); solution(n+1).pressure(1:6:end)=solution(n+1).u_g(4:6:end);
0105 
0106     %compute new thickness
0107     disp(sprintf('\n%s','   computing new thickness...'));
0108     %compute depth averaged horizontal velocity
0109     velocity_average=HorizontalVelocityDepthAverage(md,solution(n+1).u_g);
0110     fem.inputs.velocity_average=velocity_average;
0111     new_thickness=iceprognostic_core(fem.m_p,'prognostic',fem.inputs);
0112 
0113     %project collapsed thickness onto 3d mesh
0114     indx=1:6:fem.m_p.gridset.gsize; indx=indx(fem.m_p.tpart);
0115     new_thickness=project3d(md,project2d(md,new_thickness(indx),1),'node');
0116 
0117     %update surface and bed using the new thickness
0118     disp(sprintf('%s','   updating geometry...'));
0119     [new_bed,new_surface,new_thickness]=UpdateGeometry(md,new_thickness,solution(n).bed(indx),solution(n).surface(indx));
0120 
0121     %project onto 3d mesh
0122     [solution(n+1).mesh solution(n+1).u_g solution(n+1).t_g solution(n+1).pressure]=UpdateMesh(md,fem.m_dv.grids,solution(n).mesh,solution(n+1).u_g,solution(n+1).t_g,solution(n+1).pressure,new_bed,new_thickness);
0123     solution(n+1).thickness=zeros(fem.m_p.gridset.gsize,1);solution(n+1).thickness(1:6:fem.m_p.gridset.gsize)=new_thickness;
0124     solution(n+1).bed=zeros(fem.m_p.gridset.gsize,1);solution(n+1).bed(1:6:fem.m_p.gridset.gsize)=new_bed;
0125     solution(n+1).surface=zeros(fem.m_p.gridset.gsize,1);solution(n+1).surface(1:6:fem.m_p.gridset.gsize)=new_surface;
0126 
0127     %figure out if time stepping is good
0128     disp(sprintf('%s','   checking time stepping...'));
0129     [back,dt,time]=TimeStepping(md,solution,dt,time);
0130     if back,
0131         continue;
0132     end
0133 
0134     %update grids
0135     disp(sprintf('%s','   updating grid positions...'));
0136     grids_t=UpdateGridPosition(md,fem.m_t.grids,new_bed,new_thickness);
0137     grids_dhu=UpdateGridPosition(md,fem.m_dhu.grids,new_bed,new_thickness);
0138     grids_dbv=UpdateGridPosition(md,fem.m_dbv.grids,new_bed,new_thickness);
0139     grids_dv=UpdateGridPosition(md,fem.m_dv.grids,new_bed,new_thickness);
0140     grids_p=UpdateGridPosition(md,fem.m_p.grids,new_bed,new_thickness);
0141 
0142     %update time and counter
0143     time=time+dt;
0144     n=n+1;
0145 
0146 end
0147 
0148 %load results onto model
0149 indx=1:6:fem.m_dv.gridset.gsize; indx=indx(fem.m_dv.tpart);
0150 indy=2:6:fem.m_dv.gridset.gsize; indy=indy(fem.m_dv.tpart);
0151 indz=3:6:fem.m_dv.gridset.gsize; indz=indz(fem.m_dv.tpart);
0152 for i=1:length(solution),
0153     solution2(i).vx=solution(i).u_g(indx)*md.yts;
0154     solution2(i).vy=solution(i).u_g(indy)*md.yts;
0155     solution2(i).vz=solution(i).u_g(indz)*md.yts;
0156     solution2(i).vel=sqrt(solution2(i).vx.^2+solution2(i).vy.^2+solution2(i).vz.^2);
0157     solution2(i).pressure=solution(i).pressure(indx)/10^5;
0158     solution2(i).temperature=solution(i).t_g(indx);
0159     solution2(i).melting=solution(i).melting_g(indx);
0160     solution2(i).thickness=solution(i).thickness(indx);
0161     solution2(i).surface=solution(i).surface(indx);
0162     solution2(i).bed=solution(i).bed(indx);
0163     solution2(i).time=solution(i).time/md.yts;
0164 end
0165 
0166 md.transient_results=solution2;

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