Changeset 2223

Timestamp:

09/11/09 16:50:20 (16 years ago)

Author:

Eric.Larour

Message:

Added new steady state control test.
Added reset of stabilization in thermal_core_nonlinear, so that frozen penalties are unfrozen at each
step of a thermal analysis.

Location:

issm/trunk

Files:

• src/c/objects/Pengrid.cpp (modified) (4 diffs)
• src/c/parallel/control_core.cpp (modified) (2 diffs)
• src/c/parallel/thermal_core_nonlinear.cpp (modified) (1 diff)
• src/m/solutions/cielo/control_core.m (modified) (2 diffs)
• test/Verification/PigControlMethodDragP3d_22/Archive_steadystate_none_0_1_absolute_1.mat (added)
• test/Verification/PigControlMethodDragP3d_22/Pig.par (modified) (1 diff)
• test/Verification/PigControlMethodDragP3d_22/configuration.m (modified) (1 diff)
• test/Verification/PigControlMethodDragP3d_22/testpresolve.m (added)

issm/trunk/src/c/objects/Pengrid.cpp

@@ -60 +60 @@
         printf("   matpar_offset=%i\n",matpar_offset);
         
-        if(node)node->Echo();
-        if(matpar)matpar->Echo();
         return;
 }
@@ -560 +558 @@
         int  dofs1[1]={0};
         int  unstable=0;
+        int  reset_penalties=0;
         
         ParameterInputs* inputs=NULL;
@@ -579 +578 @@
         found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
         if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
-
+        
+        found=inputs->Recover("reset_penalties",&reset_penalties);
+
+        if(reset_penalties)zigzag_counter=0;
 
         //Compute pressure melting point
@@ -608 +610 @@
         /*If penalty keeps zigzagging more than 5 times: */
         if(stabilize_constraints){
-                if(zigzag_counter>5){
+                if(zigzag_counter>stabilize_constraints){
                         unstable=0;
                         active=1;

issm/trunk/src/c/parallel/control_core.cpp

@@ -94 +94 @@
                 inputs->Add("fit",fit[n]);
 
-                /*Update parameters: */
-                UpdateFromInputsx(fem_model->elements,fem_model->nodes,fem_model->loads, fem_model->materials,inputs);
-
-                _printf_("%s\n","      computing gradJ...");
-                gradjcompute_results=new DataSet(ResultsEnum()); 
-                gradjcompute_core(gradjcompute_results,model, inputs);
-                gradjcompute_results->FindResult(&grad_g,"grad_g");
-                delete gradjcompute_results;
-                _printf_("%s\n","      done.");
-
-                _printf_("%s\n","      normalizing directions...");
-                Orthx(&new_grad_g,grad_g,grad_g_old);
-                VecFree(&grad_g); VecFree(&grad_g_old); 
-                grad_g_old=new_grad_g;
-                VecToMPISerial(&grad_g_double,new_grad_g);
-                _printf_("%s\n","      done.");
-
-                _printf_("%s\n","      optimizing along gradient direction...");
-                optargs.model=model;
-                optargs.param_g=param_g; optargs.grad_g=grad_g_double; optargs.inputs=inputs;optargs.n=n;
-                optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx; optpars.maxiter=(int)maxiter[n];optpars.cmjump=cmjump[n];
-                BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
-                _printf_("%s\n","      done.");
-
-                _printf_("%s\n","      updating parameter using optimized search scalar...");
-                for(i=0;i<numberofnodes;i++)param_g[i]=param_g[i]+search_scalar*optscal[n]*grad_g_double[i];
-                _printf_("%s\n","      done.");
-
-                _printf_("%s\n","      constraining the new distribution...");    
-                ControlConstrainx(param_g,numberofnodes,mincontrolconstraint,maxcontrolconstraint,control_type);
-                _printf_("%s\n","      done.");
-
                 /*In case we are running a steady state control method, compute new temperature field using new parameter 
                  * distribution: */
@@ -133 +101 @@
                 delete steadystate_results;
                 inputs->Add("temperature",t_g,1,numberofnodes);
-
+        
+                /*Update parameters: */
+                UpdateFromInputsx(fem_model->elements,fem_model->nodes,fem_model->loads, fem_model->materials,inputs);
+
+                _printf_("%s\n","      computing gradJ...");
+                gradjcompute_results=new DataSet(ResultsEnum()); 
+                gradjcompute_core(gradjcompute_results,model, inputs);
+                gradjcompute_results->FindResult(&grad_g,"grad_g");
+                delete gradjcompute_results;
+                _printf_("%s\n","      done.");
+
+                _printf_("%s\n","      normalizing directions...");
+                Orthx(&new_grad_g,grad_g,grad_g_old);
+                VecFree(&grad_g); VecFree(&grad_g_old); 
+                grad_g_old=new_grad_g;
+                VecToMPISerial(&grad_g_double,new_grad_g);
+                _printf_("%s\n","      done.");
+
+                _printf_("%s\n","      optimizing along gradient direction...");
+                optargs.model=model;
+                optargs.param_g=param_g; optargs.grad_g=grad_g_double; optargs.inputs=inputs;optargs.n=n;
+                optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx; optpars.maxiter=(int)maxiter[n];optpars.cmjump=cmjump[n];
+                BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
+                _printf_("%s\n","      done.");
+
+                _printf_("%s\n","      updating parameter using optimized search scalar...");
+                for(i=0;i<numberofnodes;i++)param_g[i]=param_g[i]+search_scalar*optscal[n]*grad_g_double[i];
+                _printf_("%s\n","      done.");
+
+                _printf_("%s\n","      constraining the new distribution...");    
+                ControlConstrainx(param_g,numberofnodes,mincontrolconstraint,maxcontrolconstraint,control_type);
+                _printf_("%s\n","      done.");
+
+        
+                
                 _printf_("%s%i%s%g\n","      value of misfit J after optimization #",n+1,": ",J[n]);
 

issm/trunk/src/c/parallel/thermal_core_nonlinear.cpp

@@ -58 +58 @@
 
                 if(debug)_printf_("%s\n","starting direct shooting method");
+
+                if(count==1) inputs->Add("reset_penalties",1);
+                else inputs->Add("reset_penalties",0);
 
                 /*Update parameters: */

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

@@ -33 +33 @@
 
         disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(model.parameters.nsteps)]));
+
+        %In case we are running a steady state control method, compute new temperature field using new parameter distribution: 
+        if model.parameters.control_steady;
+                steadystate_results=steadystate_core(models,inputs); t_g=results.t_g; 
+                inputs=add(inputs,'temperature',t_g,'doublevec',1,model.parameters.numberofnodes);
+        end
 
         %update inputs with new fit
@@ -75 +81 @@
         param_g=ControlConstrain(param_g,model.parameters);
 
-        %In case we are running a steady state control method, compute new temperature field using new parameter distribution: 
-        if model.parameters.control_steady;
-                steadystate_results=steadystate_core(models,inputs); t_g=results.t_g; 
-                inputs=add(inputs,'temperature',t_g,'doublevec',1,model.parameters.numberofnodes);
-        end
-
+        
         disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]);
 

issm/trunk/test/Verification/PigControlMethodDragP3d_22/Pig.par

@@ -52 +52 @@
 md.maxiter=10*ones(md.nsteps,1);
 md.cmjump=0.3*ones(md.nsteps,1);
+md.debug=0;
+md.stabilize_constraints=5;

issm/trunk/test/Verification/PigControlMethodDragP3d_22/configuration.m

@@ -26 +26 @@
                                  {'diagnostic',  'none',      0 ,     1,      'relative',    1    };...
                                  {'diagnostic',  'none',      0 ,     1,      'logarithmic', 1    };...
+                                 {'steadystate',  'none',      0 ,     1,      'absolute', 1    };...
         };