Changeset 1205

Timestamp:

07/01/09 16:30:57 (16 years ago)

Author:

Mathieu Morlighem

Message:

Stokes later... (pas trop vite en besogne...)

File:

• issm/trunk/src/m/solutions/cielo/control_core.m (modified) (4 diffs)

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

@@ -8 +8 @@
 %recover models
 m_dh=models.dh;
-m_dv=models.dv;
-m_ds=models.ds;
-m_sl=models.sl;
 
 %recover parameters common to all solutions
@@ -25 +22 @@
 options=ControlOptions(m_dh.parameters);
 
-%Take car of Stokes : compute slope and get spc once for all
-if isstokes,
-
-        %bed slope
-        displaystring(debug,'\n%s',['      computing bed slope (x and y derivatives)...']);
-        slopex=diagnostic_core_linear(m_sl,inputs,'slope_compute','bedx');
-        slopey=diagnostic_core_linear(m_sl,inputs,'slope_compute','bedy');
-        slopex=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopex,'slopex',0);
-        slopey=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopey,'slopey',0);
-        inputs=add(inputs,'bedslopex',slopex,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes);
-        inputs=add(inputs,'bedslopey',slopey,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes);
-
-        %3d velocity
-        displaystring(debug,'\n%s',['      get vertical spcs: launch diagnostic...']);
-        displaystring(debug,'\n%s',['         computing horizontal velocities...']);
-        u_g=diagnostic_core_nonlinear(m_dh,inputs,'diagnostic','horiz');
-        displaystring(debug,'\n%s',['         extruding horizontal velocities...']);
-        u_g_horiz=FieldExtrude(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,u_g,'velocity',1);
-        displaystring(debug,'\n%s',['         computing vertical velocities...']);
-        inputs=add(inputs,'velocity',u_g_horiz,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes);
-        u_g_vert=diagnostic_core_linear(m_dv,inputs,'diagnostic','vert');
-        displaystring(debug,'\n%s',['         combining horizontal and vertical velocities...']);
-        u_g=zeros(m_dh.parameters.numberofnodes*3,1);
-        u_g(dofsetgen([1,2],3,m_dh.parameters.numberofnodes*3))=u_g_horiz;
-        u_g(dofsetgen([3],3,m_dh.parameters.numberofnodes*3))=u_g_vert;
-        displaystring(debug,'\n%s',['         computing pressure according to MacAyeal...']);
-        p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,m_dh.parameters,inputs);
-        p_g=p_g/m_ds.parameters.stokesreconditioning;
-        %recombine u_g and p_g: 
-        u_g_stokes=zeros(m_ds.nodesets.gsize,1);
-        u_g_stokes(dofsetgen([1,2,3],4,m_ds.nodesets.gsize))=u_g;
-        u_g_stokes(dofsetgen([4],4,m_ds.nodesets.gsize))=p_g;
-        inputs=add(inputs,'velocity',u_g_stokes,'doublevec',4,m_ds.parameters.numberofnodes);
-        displaystring(debug,'\n%s',['         update boundary conditions for stokes using velocities previously computed...']);
-        m_ds.y_g=zeros(m_ds.nodesets.gsize,1);
-        m_ds.y_g(dofsetgen([1,2,3],4,m_ds.nodesets.gsize))=u_g;
-        [m_ds.ys m_ds.ys0]=Reducevectorgtos(m_ds.y_g,m_ds.nodesets);
-        displaystring(debug,'\n%s',['         computing stokes velocities and pressure ...']);
-        u_g=diagnostic_core_nonlinear(m_ds,inputs,'diagnostic','stokes');
-        inputs=add(inputs,'velocity',u_g,'doublevec',4,m_ds.parameters.numberofnodes);
-end
-
 for n=1:m_dh.parameters.nsteps,
 
@@ -79 +34 @@
 
         %Update inputs in datasets
-        if isstokes,
-                [m_ds.elements,m_ds.nodes,m_ds.loads,m_ds.materials]=UpdateFromInputs(m_ds.elements,m_ds.nodes,m_ds.loads,m_ds.materials,inputs);
-        else
-                [m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials]=UpdateFromInputs(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,inputs);
-        end
+        [m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials]=UpdateFromInputs(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,inputs);
 
         displaystring(debug,'\n%s',['      computing gradJ...']);
-        if isstokes,
-                [u_g c(n).grad_g]=GradJCompute(m_ds,inputs,'diagnostic','stokes');
-                inputs=add(inputs,'velocity',u_g,'doublevec',4,m_ds.parameters.numberofnodes);
-        else
-                [u_g c(n).grad_g]=GradJCompute(m_dh,inputs,'diagnostic','horiz');
-                inputs=add(inputs,'velocity',u_g,'doublevec',2,m_dh.parameters.numberofnodes);
-        end
+        [u_g c(n).grad_g]=GradJCompute(m_dh,inputs,'diagnostic','horiz');
+        inputs=add(inputs,'velocity',u_g,'doublevec',2,m_dh.parameters.numberofnodes);
 
         displaystring(debug,'\n%s',['      normalizing directions...']);
@@ -107 +53 @@
 
         displaystring(debug,'\n%s',['      optimizing along gradient direction...']);
-        if isstokes,
-                [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m_dh,inputs,param_g,c(n).grad_g,n,'diagnostic','horiz');
-        else
-                [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m_ds,inputs,param_g,c(n).grad_g,n,'diagnostic','horiz');
-        end
+        [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m_dh,inputs,param_g,c(n).grad_g,n,'diagnostic','horiz');
 
         displaystring(debug,'\n%s',['      updating parameter using optimized search scalar...']);