Changeset 3203

Timestamp:

03/05/10 14:51:26 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added field weights for CM

Location:

issm/trunk/src

Files:

• c/ModelProcessorx/Control/CreateParametersControl.cpp (modified) (2 diffs)
• c/ModelProcessorx/IoModel.cpp (modified) (2 diffs)
• c/ModelProcessorx/IoModel.h (modified) (1 diff)
• c/objects/Tria.cpp (modified) (7 diffs)
• c/parallel/diagnostic.cpp (modified) (3 diffs)
• c/parallel/steadystate.cpp (modified) (2 diffs)
• m/classes/@model/model.m (modified) (2 diffs)
• m/classes/public/display/displaycontrol.m (modified) (1 diff)
• m/classes/public/ismodelselfconsistent.m (modified) (5 diffs)
• m/classes/public/marshall.m (modified) (1 diff)
• m/solutions/jpl/diagnostic.m (modified) (1 diff)

issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp

@@ -22 +22 @@
 
         double* fit=NULL;
+        double* weights=NULL;
         double* cm_jump=NULL;
         double* optscal=NULL;
@@ -176 +177 @@
                 xfree((void**)&iomodel->optscal);
                 xfree((void**)&iomodel->maxiter);
+
+                /*Now, recover weights as vectors: */
+                IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights");
+
+                count++;
+                param= new Param(count,"weights",DOUBLEVEC);
+                param->SetDoubleVec(iomodel->weights,iomodel->numberofnodes);
+                parameters->AddObject(param);
+
+                xfree((void**)&iomodel->weights);
 
                 /*Get vx, vx_obs, vy, vy_obs, and the parameter value: */

issm/trunk/src/c/ModelProcessorx/IoModel.cpp

@@ -108 +108 @@
         /*!solution parameters: */
         iomodel->fit=NULL;
+        iomodel->weights=NULL;
         iomodel->cm_jump=NULL;
         iomodel->meanvel=0;
@@ -252 +253 @@
         xfree((void**)&iomodel->n);
         xfree((void**)&iomodel->fit);
+        xfree((void**)&iomodel->weights);
         xfree((void**)&iomodel->cm_jump);
         xfree((void**)&iomodel->optscal);

issm/trunk/src/c/ModelProcessorx/IoModel.h

@@ -124 +124 @@
         /*solution parameters: */
         double* fit;
+        double* weights;
         double* cm_jump;
         int     artificial_diffusivity;

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

@@ -2858 +2858 @@
         int          numberofdofspernode;
         int          dofs2[2]={0,1};
+        int          dofs1[1]={0};
 
         /* grid data: */
@@ -2868 +2869 @@
         double dux_list[numgrids];
         double duy_list[numgrids];
+        double weights_list[numgrids];
 
         /* gaussian points: */
@@ -2921 +2923 @@
         if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){
                 throw ErrorException(__FUNCT__,"missing velocity input parameter");
+        }
+        if(!inputs->Recover("weights",&weights_list[0],1,dofs1,numgrids,(void**)nodes)){
+                throw ErrorException(__FUNCT__,"missing weights input parameter");
         }
 
@@ -3042 +3047 @@
         }
 
+        /*Apply weights to DU*/
+        for (i=0;i<numgrids;i++){
+                dux_list[i]=weights_list[i]*dux_list[i];
+                duy_list[i]=weights_list[i]*duy_list[i];
+        }
+
         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
@@ -4491 +4502 @@
         double obs_vy_list[numgrids];
         double misfit_list[numgrids];
+        double weights_list[numgrids];
 
         /* gaussian points: */
@@ -4535 +4547 @@
         if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){
                 throw ErrorException(__FUNCT__,"missing velocity input parameter");
+        }
+        if(!inputs->Recover("weights",&weights_list[0],1,dofs1,numgrids,(void**)nodes)){
+                throw ErrorException(__FUNCT__,"missing weights input parameter");
         }
 
@@ -4626 +4641 @@
         }
 
+        /*Apply weights to misfits*/
+        for (i=0;i<numgrids;i++){
+                misfit_list[i]=weights_list[i]*misfit_list[i];
+        }
+
         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);

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

@@ -41 +41 @@
         double* u_g_initial=NULL;
         double* u_g_obs=NULL;
+        double* weights=NULL;
         Param*  param=NULL;
         int      count;
@@ -100 +101 @@
         model->FindParam(&u_g_initial,NULL,NULL,"u_g",DiagnosticAnalysisEnum(),HorizAnalysisEnum());
         model->FindParam(&numberofnodes,"numberofnodes");
-        if(control_analysis | qmu_analysis){
-                model->FindParam(&u_g_obs,NULL,NULL,"u_g_obs",DiagnosticAnalysisEnum(),HorizAnalysisEnum());
-        }
 
         inputs=new ParameterInputs;
+        inputs->Add("velocity",u_g_initial,3,numberofnodes);
 
-        inputs->Add("velocity",u_g_initial,3,numberofnodes);
-        if(control_analysis | qmu_analysis){
+        if(control_analysis){
+                model->FindParam(&u_g_obs,NULL,NULL,"u_g_obs",DiagnosticAnalysisEnum(),HorizAnalysisEnum());
+                model->FindParam(&weights,NULL,NULL,"weights",DiagnosticAnalysisEnum(),HorizAnalysisEnum());
                 inputs->Add("velocity_obs",u_g_obs,2,numberofnodes);
+                inputs->Add("weights",weights,1,numberofnodes);
         }
         
@@ -174 +175 @@
         xfree((void**)&u_g_initial);
         xfree((void**)&u_g_obs);
+        xfree((void**)&weights);
         delete model;
         delete inputs;

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

@@ -44 +44 @@
         double* p_g_initial=NULL;
         double* u_g_obs=NULL;
+        double* weights=NULL;
         double  dt;
         Param*  param=NULL;
@@ -119 +120 @@
                 model->FindParam(&u_g_obs,NULL,NULL,"u_g_obs",DiagnosticAnalysisEnum(),HorizAnalysisEnum());
                 inputs->Add("velocity_obs",u_g_obs,2,numberofnodes);
+                model->FindParam(&weights,NULL,NULL,"weights",DiagnosticAnalysisEnum(),HorizAnalysisEnum());
+                inputs->Add("velocity_obs",weights,1,numberofnodes);
         }
         

issm/trunk/src/m/classes/@model/model.m

@@ -187 +187 @@
         md.control_analysis=0;
         md.control_type='';
+        md.weights=[];
         md.nsteps=0;
         md.maxiter=[];
-        md.eps_cm=0;
+        md.fit=[];
         md.tolx=0;
         md.optscal=[];
+        md.eps_cm=0;
         md.cm_noisedmp=0;
         md.cm_mindmp_value=0;
@@ -199 +201 @@
         md.cm_min=0;
         md.cm_max=0;
-        md.fit=[];
         md.cm_jump=[];
         md.cm_gradient=0;

issm/trunk/src/m/classes/public/display/displaycontrol.m

@@ -12 +12 @@
 if md.control_analysis,
         fielddisplay(md,'control_type','parameter where inverse control is carried out; ex: ''drag'', or ''B'')');
+        fielddisplay(md,'weights','weights applied to the misfit of each node');
         fielddisplay(md,'nsteps','number of optimization searches');
         fielddisplay(md,'eps_cm','misfit convergence criterion. Default is 1%, NaN if not applied');
         fielddisplay(md,'optscal','scaling factor on gradient direction during optimization, for each optimization step');
         fielddisplay(md,'fit','''absolute: 0'', ''relative: 1'', or ''logarithmic: 2''. for each optimization steps');
-        fielddisplay(md,'cm_jump','decrease threshold for misfit, default is 30%');
         fielddisplay(md,'maxiter','maximum iterations during each optimization step');
         fielddisplay(md,'tolx','minimum tolerance which will stop one optimization search');
+        fielddisplay(md,'cm_jump','decrease threshold for misfit, default is 30%');
         fielddisplay(md,'cm_noisedmp','noise dampening coefficient, 0 if not applied');
         fielddisplay(md,'cm_mindmp_value','minimum acceptable value of the inversed parameter before being dampened');

issm/trunk/src/m/classes/public/ismodelselfconsistent.m

@@ -85 +85 @@
 fields={'numberofelements','numberofgrids','x','y','z','drag','drag_type','p','q',...
         'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
-        'tolx','np','eps_res','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
+        'tolx','np','eps_res','max_nonlinear_iterations','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
 checknan(md,fields);
 
 %FIELDS >= 0 
 fields={'numberofelements','numberofgrids','elements','drag','drag_type','p','q',...
-        'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_res','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
+        'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_res','max_nonlinear_iterations','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
         'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
 checkgreater(md,fields,0);
@@ -96 +96 @@
 %FIELDS > 0
 fields={'numberofelements','numberofgrids','elements','drag_type','p',...
-        'rho_ice','rho_water','B','thickness','g','eps_res','eps_rel','eps_abs','maxiter','tolx',...
+        'rho_ice','rho_water','B','thickness','g','max_nonlinear_iterations','eps_res','eps_rel','eps_abs','maxiter','tolx',...
         'sparsity','deltaH','DeltaH','timeacc','timedec'};
 checkgreaterstrict(md,fields,0);
@@ -373 +373 @@
         checkvalues(md,{'fit'},[0 1 2 3 4]);
 
+        %WEIGHTS
+        fields={'weights'};
+        checksize(md,fields,[md.numberofgrids 1]);
+        checkgreater(md,fields,0);
+
         %OBSERVED VELOCITIES
         fields={'vx_obs','vy_obs'};
@@ -418 +423 @@
                 if isnan(fieldsize(1)),
                         if (size(eval(['md.' fields{i}]),2)~=fieldsize(2)),
-                                %LOG warning:
+                                %LOG warnings:
                                 if strcmpi(fields{i},'pressureload'),
                                         disp(' ');
@@ -433 +438 @@
                         end
                 else
+                        %LOG warnings:
+                        if strcmpi(fields{i},'weights'),
+                                disp(' ');
+                                disp('   ''weights'' is a new field that should be of length numberofgrids.');
+                                disp('      To update your model, use the following line');
+                                disp('      md.weights=ones(md.numberofgrids,1)');
+                                disp(' ');
+                        end
                         if ((size(eval(['md.' fields{i}]),1)~=fieldsize(1)) |  (size(eval(['md.' fields{i}]),2)~=fieldsize(2)))
                                 error(['model not consistent: field ' fields{i} ' size should be ' num2str(fieldsize(1)) ' x ' num2str(fieldsize(2))]);

issm/trunk/src/m/classes/public/marshall.m

@@ -100 +100 @@
 %Write solution parameters
 WriteData(fid,md.fit,'Mat','fit');
+WriteData(fid,md.weights,'Mat','weights');
 WriteData(fid,md.cm_jump,'Mat','cm_jump');
 WriteData(fid,md.yts,'Scalar','yts');

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

@@ -34 +34 @@
         if md.control_analysis,
                 inputs=add(inputs,'velocity_obs',models.dh.parameters.u_g_obs,'doublevec',2,models.dh.parameters.numberofnodes);
+                inputs=add(inputs,'weights',models.dh.parameters.weights,'doublevec',1,models.dh.parameters.numberofnodes);
         end