Changeset 727

Timestamp:

06/02/09 15:13:44 (16 years ago)

Author:

Mathieu Morlighem

Message:

corrected prognostic core (identical as parallel)

Location:

issm/trunk/src/m/solutions/cielo

Files:

• prognostic.m (modified) (3 diffs)
• prognostic_core.m (modified) (1 diff)
• thermal.m (modified) (2 diffs)

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

@@ -12 +12 @@
         md.analysis_type='prognostic'; m_p=CreateFemModel(md);
 
-        %Take only the first two dofs of m_p.parameters.u_g
-        u_g=m_p.parameters.u_g(dofsetgen([1,2],3,m_p.parameters.numberofnodes*3));
-        displaystring(md.debug,'\n%s',['depth averaging velocity...']);
-        u_g=VelocityDepthAverage(m_p.elements,m_p.nodes,m_p.loads,m_p.materials,u_g);
-
         % figure out number of dof: just for information purposes.
         md.dof=m_p.nodesets.fsize;
@@ -23 +18 @@
         displaystring(md.debug,'\n%s',['setup inputs...']);
         inputs=inputlist;
-        inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m_p.parameters.numberofnodes);
+        inputs=add(inputs,'velocity',m_p.parameters.u_g,'doublevec',3,m_p.parameters.numberofnodes);
         inputs=add(inputs,'thickness',m_p.parameters.h_g,'doublevec',1,m_p.parameters.numberofnodes);
         inputs=add(inputs,'melting',m_p.parameters.m_g,'doublevec',1,m_p.parameters.numberofnodes);
@@ -29 +24 @@
         inputs=add(inputs,'dt',m_p.parameters.dt,'double');
 
-        displaystring(md.debug,'\n%s',['computing new thickness...']);
+        displaystring(md.debug,'\n%s',['call computational core:']);
         h_g=prognostic_core(m_p,inputs,'prognostic','');
-        
-        displaystring(md.debug,'\n%s',['extrude new thickness...']);
-        h_g=ThicknessExtrude(m_p.elements,m_p.nodes,m_p.loads,m_p.materials,h_g);
-        
+
         displaystring(md.debug,'\n%s',['load results...']);
         md.results.prognostic.step=1;

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

@@ -1 @@
-function u_g=prognostic_core(m,inputs,analysis_type,sub_analysis_type)
+function h_g=prognostic_core(m,inputs,analysis_type,sub_analysis_type)
 %PROGNOSTIC_CORE - linear solution sequence
 %
 %   Usage:
-%      u_g=prognostic_core(m,inputs,analysis_type,sub_analysis_type)
+%      h_g=prognostic_core(m,inputs,analysis_type,sub_analysis_type)
 
-        %stiffness and load generation only:
-        m.parameters.kflag=1; m.parameters.pflag=1;
+        displaystring(m.parameters.debug,'\n%s',['depth averaging velocity...']);
+        %Take only the first two dofs of m.parameters.u_g
+        u_g=get(inputs,'velocity',[1 1 0 0]);
+        u_g=VelocityDepthAverage(m.elements,m.nodes,m.loads,m.materials,u_g);
+        inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m.parameters.numberofnodes);
 
-        %Update inputs in datasets
-        [m.elements,m.nodes, m.loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, m.loads,m.materials,inputs);
-        
-        %system matrices
-        [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-        [K_gg, p_g,kmax]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-        
-        %Reduce tangent matrix from g size to f size
-        [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); 
-        
-        %Reduce load from g size to f size
-        [p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets);
-        
-        %Solve  
-        [u_f]=Solver(K_ff,p_f,[],m.parameters);
-        
-        %Merge back to g set
-        [u_g]= Mergesolutionfromftog( u_f, m.Gmn, m.ys, m.nodesets ); 
+        displaystring(m.parameters.debug,'\n%s',['call computational core:']);
+        h_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type);
+
+        displaystring(m.parameters.debug,'\n%s',['extrude computed thickness on all layers:']);
+        h_g=ThicknessExtrude(m.elements,m.nodes,m.loads,m.materials,h_g);
 
 end %end function

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

@@ -36 +36 @@
                 
                 displaystring(md.debug,'\n%s',['load results...']);
-                md.results.thermal.step=1;
-                md.results.thermal.time=0;
-                md.results.thermal.temperature=t_g;
-                md.results.thermal.melting=m_g*md.yts; %from m/s to m/a
+                solution=struct('step',[],'time',[],'temperature',[],'melting',[]);
+                solution.step=1;
+                solution.time=0;
+                solution.temperature=t_g;
+                solution.melting=m_g*md.yts; %from m/s to m/a;
+                md.results.thermal=solution;
 
         else
@@ -70 +72 @@
                 
                 %Wrap up
-                solution=struct('time',{},'temperature',{},'melting',{});
+                solution=struct('step',{},'time',{},'temperature',{},'melting',{});
                 for n=1:nsteps+1,
                         solution(n).temperature=soln(n).t_g;