Changeset 968

Timestamp:

06/12/09 16:52:48 (16 years ago)

Author:

Mathieu Morlighem

Message:

extended processresults to other packages

Location:

issm/trunk/src/m/solutions/cielo

Files:

• processresults.m (modified) (2 diffs)
• prognostic.m (modified) (1 diff)
• prognostic_core.m (modified) (1 diff)
• thermal.m (modified) (1 diff)
• thermal_core.m (modified) (2 diffs)
• transient3d.m (modified) (4 diffs)

issm/trunk/src/m/solutions/cielo/processresults.m

@@ -12 +12 @@
 
 %recover models first
-if strcmpi(analysis_type,'diagnostic'),
+if (strcmpi(analysis_type,'diagnostic') | strcmpi(analysis_type,'transient')),
         m_dh=models.dh;
         m_ds=models.ds;
@@ -22 +22 @@
         ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
         isstokes=m_ds.parameters.isstokes;
-
-elseif strcmpi(analysis_type,'thermal'),
+end
+if (strcmpi(analysis_type,'thermal') | strcmpi(analysis_type,'transient')),
         m_m=models.m;
 end

issm/trunk/src/m/solutions/cielo/prognostic.m

@@ -25 +25 @@
 
         displaystring(md.debug,'\n%s',['call computational core:']);
-        rawresults=prognostic_core(models,inputs,'prognostic','');
+        results=prognostic_core(models,inputs,'prognostic','');
 
         displaystring(md.debug,'\n%s',['load results...']);
         if ~isstruct(md.results), md.results=struct(); end
-        md.results.prognostic.step=1;
-        md.results.prognostic.time=0;
-        md.results.prognostic.thickness=rawresults.h_g;
+        md.results.prognostic=processresults(results,models,'prognostic');
 
         %stop timing

issm/trunk/src/m/solutions/cielo/prognostic_core.m

@@ -7 +7 @@
         %get FE model
         m=models.p;
+        results.time=0;
+        results.step=1;
 
         displaystring(m.parameters.debug,'\n%s',['depth averaging velocity...']);

issm/trunk/src/m/solutions/cielo/thermal.m

@@ -26 +26 @@
         
         %call core
-        rawresults=thermal_core(models,inputs);
-
-        %process results
-        results=struct('step',{},'time',{},'temperature',{},'melting',{});
-        yts=models.t.parameters.yts;
-        for n=1:length(rawresults),
-                results(n).temperature=rawresults(n).t_g;
-                results(n).melting=rawresults(n).m_g*yts; %in m/year
-                results(n).time=rawresults(n).time/yts;      %in years
-                results(n).step=n;
-        end
+        results=thermal_core(models,inputs);
 
         %plug onto the model
         if ~isstruct(md.results), md.results=struct(); end
-        md.results.thermal=results;
+        md.results.thermal=processresults(results,models,'thermal');
 
         %stop timing

issm/trunk/src/m/solutions/cielo/thermal_core.m

@@ -12 +12 @@
 
         results.time=0;
+        results.step=1;
 
         displaystring(m_t.parameters.debug,'\n%s',['computing temperatures...']);
@@ -29 +30 @@
 
         %initialize temperature and melting
+        results.step=1;
+        results.time=0;
         results.t_g=t_g;
         results.m_g=m_g;
-        results.time=0;
 
         for n=1:nsteps, 
 
                 displaystring(m_t.parameters.debug,'\n%s%i/%i\n','time step: ',n,nsteps);
+                results(n+1).step=n+1;
                 results(n+1).time=n*m_t.parameters.dt;
 

issm/trunk/src/m/solutions/cielo/transient3d.m

@@ -35 +35 @@
 
 %initialize solution
-solution=struct('step',[],'time',[],'u_g',[],'p_g',[],'h_g',[],'s_g',[],'b_g',[]);
-solution.step=1;
-solution.time=0;
-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.m_g;
-solution.a_g=models.p.parameters.a_g;
+results=struct('step',[],'time',[],'u_g',[],'p_g',[],'h_g',[],'s_g',[],'b_g',[]);
+results.step=1;
+results.time=0;
+results.u_g=models.dh.parameters.u_g;
+%results.u_g=models.dh.parameters.u_g(dofsetgen([1,2],3,length(models.dh.parameters.u_g)));
+results.p_g=models.p.parameters.p_g;
+results.h_g=models.p.parameters.h_g;
+results.s_g=models.p.parameters.s_g;
+results.b_g=models.p.parameters.b_g;
+results.t_g=models.p.parameters.t_g;
+results.m_g=models.p.parameters.m_g;
+results.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,models.p.parameters.numberofnodes);
-inputs=add(inputs,'melting',solution.m_g,'doublevec',1,models.p.parameters.numberofnodes);
-inputs=add(inputs,'accumulation',solution.a_g,'doublevec',1,models.p.parameters.numberofnodes);
+inputs=add(inputs,'velocity',results.u_g,'doublevec',3,models.p.parameters.numberofnodes);
+inputs=add(inputs,'melting',results.m_g,'doublevec',1,models.p.parameters.numberofnodes);
+inputs=add(inputs,'accumulation',results.a_g,'doublevec',1,models.p.parameters.numberofnodes);
 inputs=add(inputs,'dt',models.p.parameters.dt,'double');
 
@@ -70 +70 @@
         displaystring(md.debug,'\n%s%g%s%g%s%g\n','time [yr]: ',time/yts,'    iteration number: ',n,'/',floor(finaltime/dt));
 
-        solution(n+1).step=n+1;
-        solution(n+1).time=time;
+        results(n+1).step=n+1;
+        results(n+1).time=time;
 
         %update inputs
-        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);
+        inputs=add(inputs,'thickness',results(n).h_g,'doublevec',1,models.p.parameters.numberofnodes);
+        inputs=add(inputs,'surface',results(n).s_g,'doublevec',1,models.p.parameters.numberofnodes);
+        inputs=add(inputs,'bed',results(n).b_g,'doublevec',1,models.p.parameters.numberofnodes);
+        inputs=add(inputs,'velocity',results(n).u_g,'doublevec',3,models.p.parameters.numberofnodes);
+        inputs=add(inputs,'pressure',results(n).p_g,'doublevec',1,models.p.parameters.numberofnodes);
+        inputs=add(inputs,'temperature',results(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 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);
+        [results(n+1).t_g models.t.loads melting_offset]=thermal_core_nonlinear(models.t,inputs,'thermal','transient');
+        inputs=add(inputs,'temperature',results(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(models.m,inputs,'melting','transient');
+        results(n+1).m_g=diagnostic_core_linear(models.m,inputs,'melting','transient');
 
         %Compute depth averaged temperature
-        temperature_average=FieldDepthAverage(models.t.elements,models.t.nodes,models.t.loads,models.t.materials,solution(n+1).t_g,'temperature');
+        temperature_average=FieldDepthAverage(models.t.elements,models.t.nodes,models.t.loads,models.t.materials,results(n+1).t_g,'temperature');
         inputs=add(inputs,'temperature_average',temperature_average,'doublevec',1,models.t.parameters.numberofnodes);
 
         %Deal with velocities.
         rawresults=diagnostic_core(models,inputs);
-        solution(n+1).u_g=rawresults.u_g; solution(n+1).p_g=rawresults.p_g;
+        results(n+1).u_g=rawresults.u_g; results(n+1).p_g=rawresults.p_g;
 
         %compute new thickness
         displaystring(md.debug,'\n%s',['    computing new thickness...']);
-        inputs=add(inputs,'velocity',solution(n+1).u_g,'doublevec',3,models.p.parameters.numberofnodes);
+        inputs=add(inputs,'velocity',results(n+1).u_g,'doublevec',3,models.p.parameters.numberofnodes);
         rawresults=prognostic_core(models,inputs,'prognostic','');
         new_thickness=rawresults.h_g;
@@ -106 +106 @@
         %update surface and bed using the new thickness
         displaystring(md.debug,'\n%s',['    updating geometry...']);
-        [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);
+        [new_thickness,new_bed,new_surface]=UpdateGeometry(models.p.elements,models.p.nodes,models.p.loads,models.p.materials,models.p.parameters,new_thickness,results(n).b_g,results(n).s_g);
 
-        %Record bed surface and thickness in the solution
-        solution(n+1).h_g=new_thickness;
-        solution(n+1).b_g=new_bed;
-        solution(n+1).s_g=new_surface;
+        %Record bed surface and thickness in the results
+        results(n+1).h_g=new_thickness;
+        results(n+1).b_g=new_bed;
+        results(n+1).s_g=new_surface;
 
         %figure out if time stepping is good
         %displaystring(md.debug,'   checking time stepping...'));
-        %[back,dt,time]=TimeStepping(md,solution,dt,time);
+        %[back,dt,time]=TimeStepping(md,results,dt,time);
         %if back,
         %       continue;
@@ -135 +135 @@
 end
 
-%Load results onto model
-results=struct('step',[],'time',[],'vx',[],'vy',[],'vel',[],'pressure',[],'thickness',[],'surface',[],'bed',[]);
-for i=1:length(solution),
-        results(i).step=solution(i).step;
-        results(i).time=solution(i).time/yts;
-        results(i).vx=solution(i).u_g(1:3:end)*yts;
-        results(i).vy=solution(i).u_g(2:3:end)*yts;
-        results(i).vz=solution(i).u_g(3:3:end)*yts;
-        results(i).vel=sqrt(solution(i).u_g(1:3:end).^2+solution(i).u_g(2:3:end).^2+solution(i).u_g(3:3:end).^2)*yts;
-        results(i).pressure=solution(i).p_g;
-        results(i).bed=solution(i).b_g;
-        results(i).surface=solution(i).s_g;
-        results(i).thickness=solution(i).h_g;
-        results(i).temperature=solution(i).t_g;
-        results(i).melting=solution(i).m_g;
-end
+%process results
 if ~isstruct(md.results), md.results=struct(); end
-md.results.transient=results;
+md.results.transient=processresults(results,models,'transient');