Changeset 302

Timestamp:

05/07/09 12:11:37 (16 years ago)

Author:

Eric.Larour

Message:

New framework for diagnostic solution sequence

Location:

issm/trunk/src/m/solutions/cielo

Files:

• CreateFemModel.m (modified) (1 diff)
• diagnostic.m (modified) (3 diffs)
• diagnostic_core.m (added)
• diagnostic_core_nonlinear.m (modified) (1 diff)
• loadresults.m (added)

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

@@ -8 +8 @@
 function  m=CreateFEMModel(md)
 
-        disp(sprintf('\n   reading data from model %s...',md.name));
+        displaystring(md.debug,'\n   reading data from model %s...',md.name);
         [m.elements,m.nodes,m.constraints,m.loads,m.materials,parameters]=ModelProcessor(md);
         
-        disp('   generating degrees of freedom...');
+        displaystring(md.debug,'%s','   generating degrees of freedom...');
         [m.nodes,m.part,m.tpart]=Dof(m.elements,m.nodes,parameters);
 
-        disp('   generating single point constraints...');
+        displaystring(md.debug,'%s','   generating single point constraints...');
         [m.nodes,m.yg]=SpcNodes(m.nodes,m.constraints);
 
-        disp('   generating rigid body constraints...');
+        displaystring(md.debug,'%s','   generating rigid body constraints...');
         [m.Rmg,m.nodes]=MpcNodes(m.nodes,m.constraints);
 
-        disp('   generating node sets...');
+        displaystring(md.debug,'%s','   generating node sets...');
         m.nodesets=BuildNodeSets(m.nodes);
 
-        disp('   reducing single point constraints vector...');
+        displaystring(md.debug,'%s','   reducing single point constraints vector...');
         [m.ys m.ys0]=Reducevectorgtos(m.yg,m.nodesets);
 
-        disp('   normalizing rigid body constraints matrix...');
+        displaystring(md.debug,'%s','   normalizing rigid body constraints matrix...');
         m.Gmn = NormalizeConstraints(m.Rmg,m.nodesets);
         
-        disp('   configuring element and loads...');
+        displaystring(md.debug,'%s','   configuring element and loads...');
         [m.elements,m.loads,m.nodes] = ConfigureObjects( m.elements, m.loads, m.nodes, m.materials);
 
-        disp('   processing parameters...');
+        displaystring(md.debug,'%s','   processing parameters...');
         parameters= ProcessParams(parameters,m.part);
 
-        disp('   creating parameters in matlab workspace...');
+        displaystring(md.debug,'%s','   creating parameters in matlab workspace...');
         m.parameters= ParamsToWorkspace(parameters);
 

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

@@ -1 @@
-function md=diagnostic(md)
-%DIAGNOSTIC - diagnostic solution sequence
-%
+function md=diagnostic_proto(md);
+%DIAGNOSTIC - compute the velocity field of a model
 %   Usage:
 %      md=diagnostic(md)
+%
 
         %timing
@@ -9 +9 @@
 
         %Build all models requested for diagnostic simulation
+        displaystring(md.debug,'%s',['reading diagnostic horiz model data']);
         md.analysis_type='diagnostic_horiz'; m_dh=CreateFemModel(md);
+        
+        displaystring(md.debug,'\n%s',['reading diagnostic vert model data']);
+        md.analysis_type='diagnostic_vert'; m_dv=CreateFemModel(md);
+        
+        displaystring(md.debug,'\n%s',['reading diagnostic stokes model data']);
+        md.analysis_type='diagnostic_stokes'; m_ds=CreateFemModel(md);
 
-        if strcmpi(md.type,'3d'),
-                md.analysis_type='diagnostic_vert'; m_dv=CreateFemModel(md);
-        end
-
+        displaystring(md.debug,'\n%s',['reading diagnostic hutter model data']);
+        md.analysis_type='diagnostic_hutter'; m_dhu=CreateFemModel(md);
+        
+        displaystring(md.debug,'\n%s',['reading surface and bed slope computation model data']);
+        md.analysis_type='surface_slope_compute'; m_ss=CreateFemModel(md);
+        md.analysis_type='bed_slope_compute'; m_bs=CreateFemModel(md);
+        
         % figure out number of dof: just for information purposes.
         md.dof=m_dh.nodesets.fsize; %biggest dof number
@@ -22 +32 @@
         inputs=add(inputs,'velocity',m_dh.parameters.u_g,'doublevec',3,m_dh.parameters.numberofnodes);
 
-        %Get horizontal solution. 
-        disp(sprintf('\n%s',['computing horizontal velocities...']));
+        %compute solution
+        [u_g,p_g]=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_ss,m_bs,inputs);
 
-        u_g=diagnostic_core_nonlinear(m_dh,inputs);
+        %Load results onto model
+        md=loadresults(md,u_g,p_g,m_dh,m_ds,m_dhu);
 
-        if strcmpi(md.type,'3d'),
-        
-                %extrude velocities for collapsed penta elements
-                disp(sprintf('\n%s',['extruding horizontal velocities...']));
-                u_g=VelocityExtrude(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,u_g);
-
-                disp(sprintf('\n%s',['computing vertical velocities...']));
-        
-                inputs=add(inputs,'velocity',u_g,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes);
-
-                u_g_vert=diagnostic_core_linear(m_dv,inputs);
-
-                %load results onto model: 
-                md.vx=u_g(1:2:end)*md.yts;
-                md.vy=u_g(2:2:end)*md.yts;
-                md.vz=u_g_vert*md.yts;
-                md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
-                
-                %Computation of pressure with Pattyn's assumptions (P=rho_ice*g*(s-z) in Pa)
-                md.pressure=md.rho_ice*md.g*(md.surface-md.z);
-        else
-                %load results onto model
-                md.vx=u_g(1:2:end)*md.yts;
-                md.vy=u_g(2:2:end)*md.yts;
-                md.vz=zeros(md.numberofgrids,1);
-                md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
-        
-                %Computation of pressure with Pattyn's assumptions (P=rho_ice*g*(s-z) in Pa)
-                md.pressure=md.rho_ice*md.g*(md.thickness);
-        end             
-
-        %timing
+        %stop timing
         t2=clock;
-        disp(sprintf('\n%s\n',['   solution converged in ' num2str(etime(t2,t1)) ' seconds']));
+        displaystring(md.debug,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);

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

@@ -1 @@
-function [u_g varargout]=diagnostic_core_nonlinear(m,inputs)
+function [u_g varargout]=diagnostic_core_nonlinear(m,inputs,analysis_type)
 %INPUT function [ru_g varargout]=cielodiagnostic_core_nonlinear(m,inputs)