Changeset 1850

Timestamp:

08/24/09 18:05:31 (16 years ago)

Author:

Eric.Larour

Message:

Brachning back from issm.controlstatic, by hand

Location:

issm/trunk/src/m/solutions/cielo

Files:

• control.m (modified) (2 diffs)
• gradjcompute_core.m (added)
• objectivefunctionC.m (modified) (2 diffs)
• thermalstatic_core.m (modified) (2 diffs)

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

@@ -27 +27 @@
         md.analysis_type=SlopeComputeAnalysisEnum(); md.sub_analysis_type=NoneAnalysisEnum(); models.sl=CreateFemModel(md);
 
+        displaystring(md.debug,'\n%s',['reading thermal model data']);
+        md.analysis_type=ThermalAnalysisEnum(); md.sub_analysis_type=SteadyAnalysisEnum(); models.t=CreateFemModel(md);
+
+        displaystring(md.debug,'\n%s',['reading melting model data']);
+        md.analysis_type=MeltingAnalysisEnum(); md.sub_analysis_type=SteadyAnalysisEnum(); models.m=CreateFemModel(md);
+
+
         % figure out number of dof: just for information purposes.
         md.dof=modelsize(models);
@@ -34 +41 @@
         inputs=add(inputs,'velocity',models.dh.parameters.u_g,'doublevec',3,models.dh.parameters.numberofnodes);
         inputs=add(inputs,'velocity_obs',models.dh.parameters.u_g_obs,'doublevec',2,models.dh.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*models.t.parameters.yts,'double');
 
         %compute solution

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

@@ -1 @@
-function J =objectivefunctionC(search_scalar,m,inputs,p_g,grad_g,n,analysis_type,sub_analysis_type);
-        
+function J =objectivefunctionC(search_scalar,models,inputs,p_g,grad_g,n,analysis_type,sub_analysis_type);
+
+%recover active model.
+m=models.active;
+
 %recover some parameters
 optscal=m.parameters.optscal(n);
 fit=m.parameters.fit(n);
 control_type=m.parameters.control_type;
+thermalstatic=m.parameters.thermalstatic;
 
 %Update along gradient using scalar supplied by fmincon optimization routine
@@ -13 +17 @@
 
 %Run diagnostic with updated parameters.
-u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
+if ~thermalstatic, 
+        %do a simple diagnostic, with the current temperature profile, do not look for steady state.
+        u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
+else
+        %do a full thermal mechanical steady state converged computation, much slower!
+        results=thermalstatic_core(models,inputs); u_g=results.u_g; 
+
+        %u_g ships with 3 or 4 dofs, we only want the horizontal ones!
+        if ~models.dh.parameters.isstokes, 
+                u_g=u_g(dofsetgen([1,2],3,m.parameters.numberofnodes*3));
+        end
+end
+
+%add velocity to inputs.
 inputs=add(inputs,'velocity',u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
 

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

@@ -16 +16 @@
 
 %recover parameters common to all solutions
-debug=m_dhu.parameters.debug;debug=1;
+debug=m_dhu.parameters.debug;
 dim=m_dhu.parameters.dim;
 eps_rel=m_dhu.parameters.eps_rel;
@@ -40 +40 @@
         %first compute temperature at steady state.
         if (step>1),
-                inputs=add(inputs,'velocity',results.u_g,'doublevec',ndof,m_t.parameters.numberofnodes);
+                inputs=add(inputs,'velocity',results_diagnostic.u_g,'doublevec',ndof,m_t.parameters.numberofnodes);
         end
         results_thermal=thermal_core(models,inputs);