Index: /issm/trunk/src/m/solutions/cielo/loadresults.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/loadresults.m (revision 931) +++ /issm/trunk/src/m/solutions/cielo/loadresults.m (revision 932) @@ -1,3 +1,3 @@ -function md=loadresults(md,models,results); +function md=loadresults(md,models,rawresults); %LOADRESULTS - load results onto model % @@ -6,5 +6,5 @@ % % Usage: -% md=Loadresults(md,results,m_dh,m_ds,m_dhu) +% md=Loadresults(md,models,rawresults) %recover models @@ -21,6 +21,6 @@ %recover results -u_g=results.u_g; -p_g=results.p_g; +u_g=rawresults.u_g; +p_g=rawresults.p_g; if ~isstruct(md.results), md.results=struct(); end Index: /issm/trunk/src/m/solutions/cielo/thermal.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/thermal.m (revision 931) +++ /issm/trunk/src/m/solutions/cielo/thermal.m (revision 932) @@ -10,21 +10,21 @@ %Build all models requested for diagnostic simulation displaystring(md.debug,'%s',['reading thermal model data']); - md.analysis_type='thermal'; m_t=CreateFemModel(md); + md.analysis_type='thermal'; models.t=CreateFemModel(md); displaystring(md.debug,'%s',['reading melting model data']); - md.analysis_type='melting'; m_m=CreateFemModel(md); + md.analysis_type='melting'; models.m=CreateFemModel(md); % figure out number of dof: just for information purposes. - md.dof=m_t.nodesets.fsize; + md.dof=models.t.nodesets.fsize; %initialize inputs displaystring(md.debug,'\n%s',['setup inputs...']); inputs=inputlist; - inputs=add(inputs,'velocity',m_t.parameters.u_g,'doublevec',3,m_t.parameters.numberofnodes); - inputs=add(inputs,'pressure',m_t.parameters.p_g,'doublevec',1,m_t.parameters.numberofnodes); - inputs=add(inputs,'dt',m_t.parameters.dt,'double'); + inputs=add(inputs,'velocity',models.t.parameters.u_g,'doublevec',3,models.t.parameters.numberofnodes); + inputs=add(inputs,'pressure',models.t.parameters.p_g,'doublevec',1,models.t.parameters.numberofnodes); + inputs=add(inputs,'dt',models.t.parameters.dt,'double'); %call core - results=thermal_core(m_t,m_m,inputs,md.sub_analysis_type); + results=thermal_core(models.t,models.m,inputs,md.sub_analysis_type); %plug onto the model Index: /issm/trunk/src/m/solutions/cielo/transient2d.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/transient2d.m (revision 931) +++ /issm/trunk/src/m/solutions/cielo/transient2d.m (revision 932) @@ -9,20 +9,20 @@ %Build all models requested displaystring(md.debug,'%s',['reading diagnostic horiz model data']); -md.analysis_type='diagnostic'; md.sub_analysis_type='horiz'; m_dh=CreateFemModel(md); +md.analysis_type='diagnostic'; md.sub_analysis_type='horiz'; models.dh=CreateFemModel(md); displaystring(md.debug,'\n%s',['reading diagnostic vert model data']); -md.analysis_type='diagnostic'; md.sub_analysis_type='vert'; m_dv=CreateFemModel(md); +md.analysis_type='diagnostic'; md.sub_analysis_type='vert'; models.dv=CreateFemModel(md); displaystring(md.debug,'\n%s',['reading diagnostic stokes model data']); -md.analysis_type='diagnostic'; md.sub_analysis_type='stokes'; m_ds=CreateFemModel(md); +md.analysis_type='diagnostic'; md.sub_analysis_type='stokes'; models.ds=CreateFemModel(md); displaystring(md.debug,'\n%s',['reading diagnostic hutter model data']); -md.analysis_type='diagnostic'; md.sub_analysis_type='hutter'; m_dhu=CreateFemModel(md); +md.analysis_type='diagnostic'; md.sub_analysis_type='hutter'; models.dhu=CreateFemModel(md); displaystring(md.debug,'\n%s',['reading surface and bed slope computation model data']); -md.analysis_type='slope_compute'; md.sub_analysis_type=''; m_sl=CreateFemModel(md); +md.analysis_type='slope_compute'; md.sub_analysis_type=''; models.sl=CreateFemModel(md); displaystring(md.debug,'%s',['reading prognostic model data']); -md.analysis_type='prognostic'; m_p=CreateFemModel(md); +md.analysis_type='prognostic'; models.p=CreateFemModel(md); %initialize solution @@ -30,26 +30,26 @@ solution.step=1; solution.time=0; -solution.u_g=m_dh.parameters.u_g(dofsetgen([1,2],3,length(m_dh.parameters.u_g))); +solution.u_g=models.dh.parameters.u_g(dofsetgen([1,2],3,length(models.dh.parameters.u_g))); solution.p_g=[]; -solution.h_g=m_p.parameters.h_g; -solution.s_g=m_p.parameters.s_g; -solution.b_g=m_p.parameters.b_g; +solution.h_g=models.p.parameters.h_g; +solution.s_g=models.p.parameters.s_g; +solution.b_g=models.p.parameters.b_g; %initialize inputs displaystring(md.debug,'\n%s',['setup inputs...']); inputs=inputlist; -inputs=add(inputs,'velocity',solution.u_g,'doublevec',2,m_p.parameters.numberofnodes); -inputs=add(inputs,'melting',m_p.parameters.m_g,'doublevec',1,m_p.parameters.numberofnodes); -inputs=add(inputs,'accumulation',m_p.parameters.a_g,'doublevec',1,m_p.parameters.numberofnodes); -inputs=add(inputs,'dt',m_p.parameters.dt,'double'); +inputs=add(inputs,'velocity',solution.u_g,'doublevec',2,models.p.parameters.numberofnodes); +inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.numberofnodes); +inputs=add(inputs,'accumulation',models.p.parameters.a_g,'doublevec',1,models.p.parameters.numberofnodes); +inputs=add(inputs,'dt',models.p.parameters.dt,'double'); % figure out number of dof: just for information purposes. -md.dof=modelsize(m_dh,m_dv,m_ds,m_dhu,m_sl); +md.dof=modelsize(models); %first time step is given by model. -dt=m_p.parameters.dt; -finaltime=m_p.parameters.ndt; +dt=models.p.parameters.dt; +finaltime=models.p.parameters.ndt; time=dt; -yts=m_p.parameters.yts; +yts=models.p.parameters.yts; n=1; %counter @@ -62,21 +62,22 @@ %update inputs - inputs=add(inputs,'thickness',solution(n).h_g,'doublevec',1,m_p.parameters.numberofnodes); - inputs=add(inputs,'surface',solution(n).s_g,'doublevec',1,m_p.parameters.numberofnodes); - inputs=add(inputs,'bed',solution(n).b_g,'doublevec',1,m_p.parameters.numberofnodes); + inputs=add(inputs,'thickness',solution(n).h_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'surface',solution(n).s_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'bed',solution(n).b_g,'doublevec',1,models.p.parameters.numberofnodes); %Deal with velocities. %Get horizontal solution. - [solution(n+1).u_g solution(n+1).p_g]=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_sl,inputs); + resultsdiag=diagnostic_core(models,inputs); + solution(n+1).u_g=resultsdiag.u_g; solution(n+1).p_g=resultsdiag.p_g; %compute new thickness disp(sprintf('%s',' computing new thickness...')); - inputs=add(inputs,'velocity',solution(n+1).u_g,'doublevec',2,m_p.parameters.numberofnodes); - new_thickness=prognostic_core(m_p,inputs,'prognostic',''); + inputs=add(inputs,'velocity',solution(n+1).u_g,'doublevec',2,models.p.parameters.numberofnodes); + new_thickness=prognostic_core(models.p,inputs,'prognostic',''); %update surface and bed using the new thickness disp(sprintf('%s',' updating geometry...')); - [new_thickness,new_bed,new_surface]=UpdateGeometry(m_p.elements,m_p.nodes,m_p.loads,m_p.materials,m_p.parameters,new_thickness,solution(n).b_g,solution(n).s_g); + [new_thickness,new_bed,new_surface]=UpdateGeometry(models.p.elements,models.p.nodes,models.p.loads,models.p.materials,models.p.parameters,new_thickness,solution(n).b_g,solution(n).s_g); %Record bed surface and thickness in the solution Index: /issm/trunk/src/m/solutions/cielo/transient3d.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/transient3d.m (revision 931) +++ /issm/trunk/src/m/solutions/cielo/transient3d.m (revision 932) @@ -9,27 +9,27 @@ %Build all models related to diagnostic displaystring(md.debug,'%s',['reading diagnostic horiz model data']); -md.analysis_type='diagnostic'; md.sub_analysis_type='horiz'; m_dh=CreateFemModel(md); +md.analysis_type='diagnostic'; md.sub_analysis_type='horiz'; models.dh=CreateFemModel(md); displaystring(md.debug,'\n%s',['reading diagnostic vert model data']); -md.analysis_type='diagnostic'; md.sub_analysis_type='vert'; m_dv=CreateFemModel(md); +md.analysis_type='diagnostic'; md.sub_analysis_type='vert'; models.dv=CreateFemModel(md); displaystring(md.debug,'\n%s',['reading diagnostic stokes model data']); -md.analysis_type='diagnostic'; md.sub_analysis_type='stokes'; m_ds=CreateFemModel(md); +md.analysis_type='diagnostic'; md.sub_analysis_type='stokes'; models.ds=CreateFemModel(md); displaystring(md.debug,'\n%s',['reading diagnostic hutter model data']); -md.analysis_type='diagnostic'; md.sub_analysis_type='hutter'; m_dhu=CreateFemModel(md); +md.analysis_type='diagnostic'; md.sub_analysis_type='hutter'; models.dhu=CreateFemModel(md); displaystring(md.debug,'\n%s',['reading surface and bed slope computation model data']); -md.analysis_type='slope_compute'; md.sub_analysis_type=''; m_sl=CreateFemModel(md); +md.analysis_type='slope_compute'; md.sub_analysis_type=''; models.sl=CreateFemModel(md); displaystring(md.debug,'%s',['reading prognostic model data']); -md.analysis_type='prognostic'; m_p=CreateFemModel(md); +md.analysis_type='prognostic'; models.p=CreateFemModel(md); %Build all models related to thermal displaystring(md.debug,'%s',['reading thermal model data']); -md.analysis_type='thermal'; m_t=CreateFemModel(md); +md.analysis_type='thermal'; models.t=CreateFemModel(md); displaystring(md.debug,'%s',['reading melting model data']); -md.analysis_type='melting'; m_m=CreateFemModel(md); +md.analysis_type='melting'; models.m=CreateFemModel(md); @@ -38,30 +38,30 @@ solution.step=1; solution.time=0; -solution.u_g=m_dh.parameters.u_g; -%solution.u_g=m_dh.parameters.u_g(dofsetgen([1,2],3,length(m_dh.parameters.u_g))); -solution.p_g=m_p.parameters.p_g; -solution.h_g=m_p.parameters.h_g; -solution.s_g=m_p.parameters.s_g; -solution.b_g=m_p.parameters.b_g; -solution.t_g=m_p.parameters.t_g; -solution.m_g=m_p.parameters.m_g; -solution.a_g=m_p.parameters.a_g; +solution.u_g=models.dh.parameters.u_g; +%solution.u_g=models.dh.parameters.u_g(dofsetgen([1,2],3,length(models.dh.parameters.u_g))); +solution.p_g=models.p.parameters.p_g; +solution.h_g=models.p.parameters.h_g; +solution.s_g=models.p.parameters.s_g; +solution.b_g=models.p.parameters.b_g; +solution.t_g=models.p.parameters.t_g; +solution.m_g=models.p.parameters.models.g; +solution.a_g=models.p.parameters.a_g; %initialize inputs displaystring(md.debug,'\n%s',['setup inputs...']); inputs=inputlist; -inputs=add(inputs,'velocity',solution.u_g,'doublevec',3,m_p.parameters.numberofnodes); -inputs=add(inputs,'melting',solution.m_g,'doublevec',1,m_p.parameters.numberofnodes); -inputs=add(inputs,'accumulation',solution.a_g,'doublevec',1,m_p.parameters.numberofnodes); -inputs=add(inputs,'dt',m_p.parameters.dt,'double'); +inputs=add(inputs,'velocity',solution.u_g,'doublevec',3,models.p.parameters.numberofnodes); +inputs=add(inputs,'melting',solution_g,'doublevec',1,models.p.parameters.numberofnodes); +inputs=add(inputs,'accumulation',solution.a_g,'doublevec',1,models.p.parameters.numberofnodes); +inputs=add(inputs,'dt',models.p.parameters.dt,'double'); % figure out number of dof: just for information purposes. -md.dof=modelsize(m_dh,m_dv,m_ds,m_dhu,m_sl); +md.dof=modelsize(models); %first time step is given by model. -dt=m_p.parameters.dt; -finaltime=m_p.parameters.ndt; +dt=models.p.parameters.dt; +finaltime=models.p.parameters.ndt; time=dt; -yts=m_p.parameters.yts; +yts=models.p.parameters.yts; n=1; %counter @@ -74,35 +74,36 @@ %update inputs - inputs=add(inputs,'thickness',solution(n).h_g,'doublevec',1,m_p.parameters.numberofnodes); - inputs=add(inputs,'surface',solution(n).s_g,'doublevec',1,m_p.parameters.numberofnodes); - inputs=add(inputs,'bed',solution(n).b_g,'doublevec',1,m_p.parameters.numberofnodes); - inputs=add(inputs,'velocity',solution(n).u_g,'doublevec',3,m_p.parameters.numberofnodes); - inputs=add(inputs,'pressure',solution(n).p_g,'doublevec',1,m_p.parameters.numberofnodes); - inputs=add(inputs,'temperature',solution(n).t_g,'doublevec',1,m_t.parameters.numberofnodes); + inputs=add(inputs,'thickness',solution(n).h_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'surface',solution(n).s_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'bed',solution(n).b_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'velocity',solution(n).u_g,'doublevec',3,models.p.parameters.numberofnodes); + inputs=add(inputs,'pressure',solution(n).p_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'temperature',solution(n).t_g,'doublevec',1,models.t.parameters.numberofnodes); %Deal with temperature first displaystring(md.debug,'\n%s',[' computing temperatures...']); - [solution(n+1).t_g m_t.loads melting_offset]=thermal_core_nonlinear(m_t,inputs,'thermal','transient'); - inputs=add(inputs,'temperature',solution(n+1).t_g,'doublevec',1,m_t.parameters.numberofnodes); + [solution(n+1).t_g models.t.loads melting_offset]=thermal_core_nonlinear(models.t,inputs,'thermal','transient'); + inputs=add(inputs,'temperature',solution(n+1).t_g,'doublevec',1,models.t.parameters.numberofnodes); displaystring(md.debug,'\n%s',[' computing melting...']); inputs=add(inputs,'melting_offset',melting_offset,'double'); - solution(n+1).m_g=diagnostic_core_linear(m_m,inputs,'melting','transient'); + solution(n+1).m_g=diagnostic_core_linear(models.m,inputs,'melting','transient'); %Compute depth averaged temperature - temperature_average=FieldDepthAverage(m_t.elements,m_t.nodes,m_t.loads,m_t.materials,solution(n+1).t_g,'temperature'); - inputs=add(inputs,'temperature_average',temperature_average,'doublevec',1,m_t.parameters.numberofnodes); + temperature_average=FieldDepthAverage(models.t.elements,models.t.nodes,models.t.loads,models.t.materials,solution(n+1).t_g,'temperature'); + inputs=add(inputs,'temperature_average',temperature_average,'doublevec',1,models.t.parameters.numberofnodes); %Deal with velocities. - [solution(n+1).u_g,solution(n+1).p_g]=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_sl,inputs); + resultsdiag=diagnostic_core(models,inputs); + solution(n+1).u_g=resultsdiag.u_g; solution(n+1).p_g=resultsdiag.p_g; %compute new thickness displaystring(md.debug,'\n%s',[' computing new thickness...']); - inputs=add(inputs,'velocity',solution(n+1).u_g,'doublevec',3,m_p.parameters.numberofnodes); - new_thickness=prognostic_core(m_p,inputs,'prognostic',''); + inputs=add(inputs,'velocity',solution(n+1).u_g,'doublevec',3,models.p.parameters.numberofnodes); + new_thickness=prognostic_core(models.p,inputs,'prognostic',''); %update surface and bed using the new thickness displaystring(md.debug,'\n%s',[' updating geometry...']); - [new_thickness,new_bed,new_surface]=UpdateGeometry(m_p.elements,m_p.nodes,m_p.loads,m_p.materials,m_p.parameters,new_thickness,solution(n).b_g,solution(n).s_g); + [new_thickness,new_bed,new_surface]=UpdateGeometry(models.p.elements,models.p.nodes,models.p.loads,models.p.materials,models.p.parameters,new_thickness,solution(n).b_g,solution(n).s_g); %Record bed surface and thickness in the solution @@ -120,11 +121,11 @@ %update node positions displaystring(md.debug,'\n%s',[' updating node positions...']); - nodes_dh=UpdateNodePositions(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,new_bed,new_thickness); - nodes_dv=UpdateNodePositions(m_dv.elements,m_dv.nodes,m_dv.loads,m_dv.materials,new_bed,new_thickness); - nodes_ds=UpdateNodePositions(m_ds.elements,m_ds.nodes,m_ds.loads,m_ds.materials,new_bed,new_thickness); - nodes_sl=UpdateNodePositions(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,new_bed,new_thickness); - nodes_p=UpdateNodePositions(m_p.elements,m_p.nodes,m_p.loads,m_p.materials,new_bed,new_thickness); - nodes_t=UpdateNodePositions(m_t.elements,m_t.nodes,m_t.loads,m_t.materials,new_bed,new_thickness); - nodes_m=UpdateNodePositions(m_m.elements,m_m.nodes,m_m.loads,m_m.materials,new_bed,new_thickness); + nodes_dh=UpdateNodePositions(models.dh.elements,models.dh.nodes,models.dh.loads,models.dh.materials,new_bed,new_thickness); + nodes_dv=UpdateNodePositions(models.dv.elements,models.dv.nodes,models.dv.loads,models.dv.materials,new_bed,new_thickness); + nodes_ds=UpdateNodePositions(models.ds.elements,models.ds.nodes,models.ds.loads,models.ds.materials,new_bed,new_thickness); + nodes_sl=UpdateNodePositions(models.sl.elements,models.sl.nodes,models.sl.loads,models.sl.materials,new_bed,new_thickness); + nodes_p=UpdateNodePositions(models.p.elements,models.p.nodes,models.p.loads,models.p.materials,new_bed,new_thickness); + nodes_t=UpdateNodePositions(models.t.elements,models.t.nodes,models.t.loads,models.t.materials,new_bed,new_thickness); + nodes_m=UpdateNodePositions(models.m.elements,models.m.nodes,models.m.loads,models.m.materials,new_bed,new_thickness); %update time and counter