cielotransient

PURPOSE ^

transient solution

SYNOPSIS ^

function md=cielotransient(md)

DESCRIPTION ^

transient solution

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function md=cielotransient(md)
0002 %transient solution
0003 
0004     analysis_t='thermaltransient';
0005     analysis_m='melting';
0006     analysis_dh='diagnostic_horiz';
0007     analysis_dbv='diagnostic_basevert';
0008     analysis_dv='diagnostic_vert';
0009     analysis_p='prognostic';
0010 
0011     %configure, and initialise data model and parameter defaults:
0012     SetParameterDefaults;
0013 
0014     %Initialize model velocity at 0 before model processing
0015     md.vx=zeros(md.numberofgrids,1); md.vy=zeros(md.numberofgrids,1); md.vz=zeros(md.numberofgrids,1);
0016 
0017     %read input data for all analyses
0018     disp(sprintf('\n   reading data from model %s...',md.name));
0019 
0020     %Build a model for the diagnostic simulation
0021     m_t=CreateFEMModel(md,analysis_t);
0022     m_m=CreateFEMModel(md,analysis_m);
0023     m_dh=CreateFEMModel(md,analysis_dh);
0024     m_dbv=CreateFEMModel(md,analysis_dbv);
0025     m_dv=CreateFEMModel(md,analysis_dv);
0026     m_p=CreateFEMModel(md,analysis_p);
0027 
0028     % buid partitioning vectors. Same partitioning for all analyses.
0029     tpart=m_t.tpart;
0030     part=m_t.part;
0031     indx=1:6:m_t.uset.gsize; indx=indx(tpart);
0032     indy=2:6:m_t.uset.gsize; indy=indy(tpart);
0033     indz=3:6:m_t.uset.gsize; indz=indz(tpart);
0034     indp=4:6:m_t.uset.gsize; indp=indp(tpart);
0035 
0036     % figure out number of dof: just for information purposes.
0037     md.dof=m_dh.uset.fsize; %biggest dof number
0038 
0039     nsteps=md.ndt/md.dt;
0040 
0041     %initialize (velocity,pressure) and (melting,temperature)
0042     solution.u_g=zeros(uset.gsize,1);
0043     if ~isempty(md.vx), solution.u_g(1:6:uset.gsize)=md.vx(part)/md.yts; end
0044     if ~isempty(md.vy), solution.u_g(2:6:uset.gsize)=md.vy(part)/md.yts; end
0045     if ~isempty(md.vz), solution.u_g(3:6:uset.gsize)=md.vz(part)/md.yts; end
0046     solution.pressure=zeros(uset.gsize,1); solution.pressure(1:6:uset.gsize)=md.rho_ice*md.g*(md.surface(part)-md.z(part)); %lithostatic pressure to spin up
0047 
0048     solution.t_g=zeros(uset.gsize,1);
0049     if ~isempty(md.temperature), 
0050         solution.t_g(1:6:uset.gsize)=md.temperature(part);
0051     else
0052         solution.t_g(1:6:uset.gsize)=md.meltingpoint;
0053     end
0054     if ~isempty(md.melting), solution.melting_g(1:6:uset.gsize)=md.melting(part);end;
0055 
0056     %initialize geometry and time
0057     solution.thickness=zeros(uset.gsize,1);solution.thickness(1:6:uset.gsize)=md.thickness(part);
0058     solution.surface=zeros(uset.gsize,1);solution.surface(1:6:uset.gsize)=md.surface(part);
0059     solution.bed=zeros(uset.gsize,1);solution.bed(1:6:uset.gsize)=md.bed(part);
0060     solution.time=0;
0061 
0062     %some thermal parameters
0063     thermalparams.beta=md.beta;
0064     thermalparams.meltingpoint=md.meltingpoint;
0065     thermalparams.latentheat=md.latentheat;
0066     thermalparams.heatcapacity=md.heatcapacity;
0067     thermalparams.penalty=md.penalty_thermal;
0068 
0069     for n=1:nsteps,
0070 
0071         disp(sprintf('\n%s%i/%i\n','time step: ',n,md.ndt/md.dt));
0072         
0073         solution(n+1).time=(n)*md.dt;
0074 
0075         %Deal with temperature first
0076         disp(sprintf('%s','   computing temperature...'));
0077 
0078         %Call core thermal computation
0079         inputs=struct('velocity',solution(n).u_g,'pressure',solution(n).pressure,'temperature',solution(n).t_g,'dt',md.dt,...
0080                     'thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface);
0081 
0082         [solution(n+1).t_g,m_t]=cielothermal_core(m_t,m_t.params,thermalparams,inputs,analysis_t); 
0083         
0084         %Call core melting computation
0085         disp(sprintf('%s','   computing melting...'));
0086         inputs=struct('pressure',solution(n).pressure,'temperature',solution(n+1).t_g,'dt',md.dt,...
0087         'thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface);
0088     
0089         solution(n+1).melting_g=cielomelting_core(m_m,m_m.params,thermalparams,inputs,analysis_m); 
0090         
0091         %Compute depth averaged temperature for 2d type elements.
0092         temperature_workspace=solution(n+1).t_g(indx);
0093         temperature_average2d_workspace=DepthAverage(md,temperature_workspace);
0094         temperature_average_workspace=project3d(md,temperature_average2d_workspace,'node');
0095         temperature_average=zeros(uset.gsize,1);temperature_average(1:6:uset.gsize)=temperature_average_workspace(part);
0096         
0097         if md.debug,
0098             t_g_workspace=SwitchPartitioning(solution(n+1).t_g,'workspace',part,tpart,[1]); 
0099             melting_g_workspace=SwitchPartitioning(solution(n+1).melting_g,'workspace',part,tpart,[1]); 
0100             plot(md,'data',t_g_workspace(1:6:uset.gsize),'title','Temperature (K)','view',3,'data',melting_g_workspace(1:6:uset.gsize)*md.yts,'title','Melting (m/a)','view',3);pause(2);
0101         end
0102 
0103         %Deal with velocities.
0104         disp(sprintf('%s','   computing horizontal velocities...'));
0105 
0106         m_dh.rifts=md.rifts;
0107         inputs=struct('thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface,...
0108                     'temperature',solution(n+1).t_g,'temperature_average',temperature_average);
0109         solution(n+1).u_g=cielodiagnostic_core_nonlinear(m_dh,m_dh.params,inputs,analysis_dh); 
0110 
0111         %Extrude velocities for collapsed penta elements
0112         u_g_workspace=SwitchPartitioning(solution(n+1).u_g,'workspace',part,tpart,[1 2]); u_g_workspace=CieloVelocityExtrude(md,u_g_workspace);
0113         
0114         %Plug velocities back into solution, with cluster partitioning:
0115         solution(n+1).u_g=SwitchPartitioning(u_g_workspace,'cluster',part,tpart,[1 2]);
0116         
0117         if md.debug,
0118             plot(md,'data',sqrt(u_g_workspace(1:6:uset.gsize).^2+u_g_workspace(2:6:uset.gsize))*md.yts,'view',3,'title','Horizontal velocity (m/a)');pause(2);
0119         end
0120 
0121         %compute base vertical velocities
0122         disp(sprintf('%s','   computing vertical base velocity...'));
0123         inputs=struct('velocity',solution(n+1).u_g,'thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface);
0124         u_g_basevert=cielodiagnostic_core_linear(m_dbv,m_dbv.params,inputs,analysis_dbv); 
0125 
0126         if md.debug,
0127             u_g_basevert_workspace=SwitchPartitioning(u_g_basevert,'workspace',part,tpart,[3]); 
0128             plot(md,'data',u_g_basevert_workspace(3:6:uset.gsize)*md.yts,'layer',1,'title','Basal Vertical velocity (m/a)');pause(2);
0129         end
0130     
0131         %use solution(n+1).u_g_basevert to constrain vertical velocity
0132         disp(sprintf('%s','   computing vertical velocity...'));
0133         m_dv.y_g=u_g_basevert;
0134     
0135         %reduce constraining vector to s-set (the one we solve on)
0136         m_dv.y_s = Reducevectorg( m_dv.y_g,m_dv.uset);
0137 
0138         %compute depth averaged horizontal velocity
0139         velocity_average_workspace=CieloHorizontalVelocityDepthAverage(md,u_g_workspace); velocity_average=SwitchPartitioning(velocity_average_workspace,'cluster',part,tpart,[1 2]);
0140 
0141         %compute vertical velocities
0142         inputs=struct('velocity',solution(n+1).u_g,'velocity_average',velocity_average,'thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface);
0143         u_g_vert=cielodiagnostic_core_linear(m_dv,m_dv.params,inputs,analysis_dv); 
0144         
0145         
0146         %add contribution to vertical velocity
0147         solution(n+1).u_g(3:6:uset.gsize)=u_g_vert(3:6:uset.gsize);
0148     
0149         if md.debug,
0150             u_g_workspace=SwitchPartitioning(solution(n+1).u_g,'workspace',part,tpart,[3]);
0151             plot(md,'data',u_g_workspace(3:6:uset.gsize)*md.yts,'view',3,'title','Vertical velocity (m/a)');pause(2);
0152         end
0153 
0154         %compute pressure (for now, lithostatic)
0155         disp(sprintf('%s','   computing pressure...'));
0156         solution(n+1).pressure=zeros(uset.gsize,1); solution(n+1).pressure(1:6:uset.gsize)=md.rho_ice*md.g*(md.surface(part)-md.z(part));
0157         
0158 
0159         %move surface and bed velocities to first lower layer. This will be needed by the prognostic model
0160         u_g_vert_workspace=SwitchPartitioning(u_g_vert,'workspace',part,tpart,[3]);
0161         ws_workspace=ShiftLayers(md,u_g_vert_workspace(3:6:uset.gsize),md.numlayers,1);
0162         wb_workspace=ShiftLayers(md,u_g_vert_workspace(3:6:uset.gsize),1,1);
0163         
0164         melting_g_workspace=SwitchPartitioning(solution(n+1).melting_g,'workspace',part,tpart,[1]);
0165         melting_workspace=ShiftLayers(md,melting_g_workspace(1:6:uset.gsize),1,1);
0166         
0167         accumulation_workspace=ShiftLayers(md,md.accumulation,md.numlayers,1);
0168         
0169 
0170         %plug into g set
0171         ws_g=zeros(uset.gsize,1);ws_g(1:6:uset.gsize)=ws_workspace(part);
0172         wb_g=zeros(uset.gsize,1);wb_g(1:6:uset.gsize)=wb_workspace(part);
0173         melting_g=zeros(uset.gsize,1);melting_g(1:6:uset.gsize)=melting_workspace(part);
0174         accumulation_g=zeros(uset.gsize,1);accumulation_g(1:6:uset.gsize)=accumulation_workspace(part)/md.yts;
0175     
0176         
0177         if md.debug,
0178             plot(md,'data',ws_workspace*md.yts,'title','Surface velocity projected on layer 1','view',3,'data',wb_workspace*md.yts,'title','Basal velocity projected on layer 1','view',3,...
0179             'data',melting_workspace*md.yts,'title','Bed melting projected on layer 1 (m/yr)','view',3,'data',accumulation_workspace,'title','Surface accumulation projected on layer 1 (m/yr)','view',3);pause(2);
0180         end
0181         
0182         %compute new thickness
0183         disp(sprintf('%s','   computing new thickness...'));
0184         inputs=struct('thickness',solution(n).thickness,'melting',melting_g,'accumulation',accumulation_g,...
0185                        'surface_vertical_velocity',ws_g,'basal_vertical_velocity',wb_g,'dt',md.dt,'velocity_average',velocity_average);
0186         new_thickness=cieloprognostic_core(m_p,m_p.params,inputs,analysis_p);
0187     
0188         %project collapsed thickness onto 3d mesh
0189         new_thickness_workspace=SwitchPartitioning(new_thickness,'workspace',part,tpart,[1]); 
0190         new_thickness_workspace=project3d(md,project2d(md,new_thickness_workspace(1:6:end),1),'node');
0191 
0192         
0193         %update surface and bed using the new thickness
0194         [new_bed_workspace,new_surface_workspace,]=UpdateGeometry(md,new_thickness_workspace,solution(n).bed(indx),solution(n).surface(indx));
0195 
0196         if md.debug,
0197             plot(md,'data',new_thickness_workspace,'view',3,'title','New thickness (m)','data',new_surface_workspace,'view',3,'title','New surface (m)','data',new_bed_workspace,'view',3,'title','New bed (m)');pause(2);
0198         end
0199 
0200         %project onto g set
0201         solution(n+1).thickness=zeros(uset.gsize,1);solution(n+1).thickness(1:6:uset.gsize)=new_thickness_workspace(part);
0202         solution(n+1).bed=zeros(uset.gsize,1);solution(n+1).bed(1:6:uset.gsize)=new_bed_workspace(part);
0203         solution(n+1).surface=zeros(uset.gsize,1);solution(n+1).surface(1:6:uset.gsize)=new_surface_workspace(part);
0204 
0205         %update grid positions
0206         [m_t.bgpdt,m_t.bgpdtb]=CieloUpdateGridPosition(m_t.bgpdt,m_t.bgpdtb,solution(n+1).bed,solution(n+1).thickness);
0207 
0208         [m_dh.bgpdt,m_dh.bgpdtb]=CieloUpdateGridPosition(m_dh.bgpdt,m_dh.bgpdtb,solution(n+1).bed,solution(n+1).thickness);
0209 
0210         [m_dbv.bgpdt,m_dbv.bgpdtb]=CieloUpdateGridPosition(m_dbv.bgpdt,m_dbv.bgpdtb,solution(n+1).bed,solution(n+1).thickness);
0211 
0212         [m_dv.bgpdt,m_dv.bgpdtb]=CieloUpdateGridPosition(m_dv.bgpdt,m_dv.bgpdtb,solution(n+1).bed,solution(n+1).thickness);
0213 
0214         [m_p.bgpdt,m_p.bgpdtb]=CieloUpdateGridPosition(m_p.bgpdt,m_p.bgpdtb,solution(n+1).bed,solution(n+1).thickness);
0215     end
0216 
0217 
0218 
0219 %load results onto model
0220 for i=1:length(solution),
0221     solution2(i).vx=solution(i).u_g(indx)*md.yts;
0222     solution2(i).vy=solution(i).u_g(indy)*md.yts;
0223     solution2(i).vz=solution(i).u_g(indz)*md.yts;
0224     solution2(i).vel=sqrt(solution2(i).vx.^2+solution2(i).vy.^2+solution2(i).vz.^2);
0225     solution2(i).pressure=solution(i).pressure(indx)/10^5;
0226     solution2(i).temperature=solution(i).t_g(indx);
0227     solution2(i).melting=solution(i).melting_g(indx);
0228     solution2(i).thickness=solution(i).thickness(indx);
0229     solution2(i).surface=solution(i).surface(indx);
0230     solution2(i).bed=solution(i).bed(indx);
0231     solution2(i).time=solution(i).time/md.yts;
0232 end
0233 md.transient_results=solution2;

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