Changeset 731

Timestamp:

06/02/09 15:42:05 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added thermal core

Location:

issm/trunk/src/m/solutions/cielo

Files:

• thermal.m (modified) (1 diff)
• thermal_core.m (modified) (1 diff)
• thermal_core_nonlinear.m (added)

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

@@ -25 +25 @@
         inputs=add(inputs,'dt',m_t.parameters.dt,'double');
         
-        if strcmpi(md.sub_analysis_type,'steady'),
-        
-                displaystring(md.debug,'\n%s',['computing temperatures...']);
-                [t_g m_t.loads melting_offset]=thermal_core(m_t,inputs,'thermal','steady');
-                
-                displaystring(md.debug,'\n%s',['computing melting...']);
-                inputs=add(inputs,'melting_offset',melting_offset,'double');
-                inputs=add(inputs,'temperature',t_g,'doublevec',1,m_t.parameters.numberofnodes);
-                m_g=diagnostic_core_linear(m_m,inputs,'melting','steady');
-                
-                displaystring(md.debug,'\n%s',['load results...']);
-                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;
+        %call core
+        results=thermal_core(m_t,m_m,inputs,md.sub_analysis_type);
 
-        else
-
-                %initialize temperature and melting
-                t_g=m_t.parameters.t_g;
-                m_g=m_m.parameters.m_g;
-                nsteps=m_t.parameters.ndt/m_t.parameters.dt;
-
-                %initialize temperature and melting
-                soln.t_g=t_g;
-                soln.m_g=m_g;
-                soln.time=0;
-
-                for n=1:nsteps, 
-
-                        displaystring(md.debug,'\n%s%i/%i\n','time step: ',n,nsteps);
-                        soln(n+1).time=n*m_t.parameters.dt;
-                        
-                        displaystring(md.debug,'\n%s',['    computing temperatures...']);
-                        inputs=add(inputs,'temperature',soln(n).t_g,'doublevec',1,m_t.parameters.numberofnodes);
-                        [soln(n+1).t_g m_t.loads melting_offset]=thermal_core(m_t,inputs,'thermal','transient');
-                        
-                        disp('   computing melting...');
-                        inputs=add(inputs,'temperature',soln(n).t_g,'doublevec',1,m_t.parameters.numberofnodes);
-                        inputs=add(inputs,'melting_offset',melting_offset,'double');
-                        soln(n+1).m_g=diagnostic_core_linear(m_m,inputs,'melting','transient');
-                        
-                end
-                
-                %Wrap up
-                solution=struct('step',{},'time',{},'temperature',{},'melting',{});
-                for n=1:nsteps+1,
-                        solution(n).temperature=soln(n).t_g;
-                        solution(n).melting=soln(n).m_g*md.yts; %in m/year
-                        solution(n).time=soln(n).time/md.yts;         %in years
-                        solution(n).step=n;
-                end
-                md.results.thermal=solution;
-        end
+        %plug onto the model
+        md.results.thermal=results;
 
         %stop timing

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

@@ -1 @@
-function [t_g ,loads, melting_offset]=thermal_core(m,inputs,analysis_type,sub_analysis_type)
-%THERMAL_CORE - core of thermal solution sequence.
-%   model is return together with temperature
-%
+function solution=thermal_core(m_t,m_m,inputs,sub_analysis_type)
+%THERMAL_CORE - core of thermal solution
 %
 %   Usage:
-%      [t_g ,loads, melting_offset]=thermal_core(m,inputs,analysis_type,sub_analysis_type);
-%    
-%      
+%      solution=thermal_core(m_t,m_m,inputs,sub_analysis_type)
 
-        count=1;
-        converged=0;
 
-%   we are going to return the loads, make them a variable of this routine
-        loads=m.loads;
+if strcmpi(sub_analysis_type,'steady'),
 
-        %stiffness and load generation only:
-        m.parameters.kflag=1; m.parameters.pflag=1;
+        displaystring(m_t.parameters.debug,'\n%s',['computing temperatures...']);
+        [t_g m_t.loads melting_offset]=thermal_core_nonlinear(m_t,inputs,'thermal','steady');
 
-        displaystring(m.parameters.debug,'\n%s',['   starting direct shooting method']);
-        
-        while(~converged),
+        displaystring(m_t.parameters.debug,'\n%s',['computing melting...']);
+        inputs=add(inputs,'melting_offset',melting_offset,'double');
+        inputs=add(inputs,'temperature',t_g,'doublevec',1,m_t.parameters.numberofnodes);
+        m_g=diagnostic_core_linear(m_m,inputs,'melting','steady');
 
-                %Update inputs in datasets
-                [m.elements,m.nodes, loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,inputs);
+        displaystring(m_t.parameters.debug,'\n%s',['load results...']);
+        solution=struct('step',[],'time',[],'temperature',[],'melting',[]);
+        solution.step=1;
+        solution.time=0;
+        solution.temperature=t_g;
+        solution.melting=m_g*m_t.parameters.yts; %from m/s to m/a;
 
-                %system matrices 
-                if ~m.parameters.lowmem 
-                        if count==1
-                                displaystring(m.parameters.debug,'%s',['   system matrices']);
-                                [K_gg_nopenalty, p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-                        end
-                        displaystring(m.parameters.debug,'%s',['   penalty system matrices']);
-                        [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-                else
-                        displaystring(m.parameters.debug,'%s',['   system matrices']);
-                        [K_gg , p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-                        displaystring(m.parameters.debug,'%s',['   penalty system matrices']);
-                        [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-                end
+else
 
-                %Reduce tangent matrix from g size to f size
-                [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); 
+        %initialize temperature and melting
+        t_g=m_t.parameters.t_g;
+        m_g=m_m.parameters.m_g;
+        nsteps=m_t.parameters.ndt/m_t.parameters.dt;
 
-                %Reduce load from g size to f size
-                [p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets);
+        %initialize temperature and melting
+        soln.t_g=t_g;
+        soln.m_g=m_g;
+        soln.time=0;
 
-                %Solve  
-                displaystring(m.parameters.debug,'%s%g','   condition number of stiffness matrix: ',condest(K_ff));
-                [t_f]=Solver(K_ff,p_f,[],m.parameters);
+        for n=1:nsteps, 
 
-                %Merge back to g set
-                displaystring(m.parameters.debug,'%s',['   merging solution back to g set']);
-                [t_g]= Mergesolutionfromftog( t_f, m.Gmn, m.ys, m.nodesets ); 
+                displaystring(m_t.parameters.debug,'\n%s%i/%i\n','time step: ',n,nsteps);
+                soln(n+1).time=n*m_t.parameters.dt;
 
-                %Update inputs in datasets
-                inputs=add(inputs,'temperature',t_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
-                displaystring(m.parameters.debug,'%s',['   update inputs']);
-                [m.elements,m.nodes, loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,inputs);
-        
-                %penalty constraints
-                displaystring(m.parameters.debug,'%s',['   penalty constraints']);
-                [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-        
-                if ~converged,
-                        displaystring(m.parameters.debug,'%s%i','   #unstable constraints ',num_unstable_constraints);
-                        
-                        if num_unstable_constraints<=m.parameters.min_thermal_constraints,
-                                converged=1;
-                        end
-                end
+                displaystring(m_t.parameters.debug,'\n%s',['    computing temperatures...']);
+                inputs=add(inputs,'temperature',soln(n).t_g,'doublevec',1,m_t.parameters.numberofnodes);
+                [soln(n+1).t_g m_t.loads melting_offset]=thermal_core_nonlinear(m_t,inputs,'thermal','transient');
 
-                count=count+1;
+                disp('   computing melting...');
+                inputs=add(inputs,'temperature',soln(n).t_g,'doublevec',1,m_t.parameters.numberofnodes);
+                inputs=add(inputs,'melting_offset',melting_offset,'double');
+                soln(n+1).m_g=diagnostic_core_linear(m_m,inputs,'melting','transient');
+
         end
 
+        %Wrap up
+        solution=struct('step',{},'time',{},'temperature',{},'melting',{});
+        for n=1:nsteps+1,
+                solution(n).temperature=soln(n).t_g;
+                solution(n).melting=soln(n).m_g*m_t.parameters.yts; %in m/year
+                solution(n).time=soln(n).time/m_t.parameters.yts;         %in years
+                solution(n).step=n;
+        end
 end