source: issm/oecreview/Archive/21724-22754/ISSM-22442-22443.diff@ 22755

Last change on this file since 22755 was 22755, checked in by Mathieu Morlighem, 7 years ago

CHG: added 21724-22754
File size: 14.2 KB
RevLine 
[22755]1Index: ../trunk-jpl/src/m/classes/adm1qn3inversion.m
2===================================================================
3--- ../trunk-jpl/src/m/classes/adm1qn3inversion.m (nonexistent)
4+++ ../trunk-jpl/src/m/classes/adm1qn3inversion.m (revision 22443)
5@@ -0,0 +1,193 @@
6+%ADM1QN3INVERSION class definition
7+%
8+% Usage:
9+% adm1qn3inversion=adm1qn3inversion();
10+
11+classdef adm1qn3inversion
12+ properties (SetAccess=public)
13+ iscontrol = 0
14+ control_parameters = NaN
15+ control_scaling_factors = NaN
16+ maxsteps = 0
17+ maxiter = 0
18+ dxmin = 0
19+ gttol = 0
20+ cost_functions = NaN
21+ cost_functions_coefficients = NaN
22+ min_parameters = NaN
23+ max_parameters = NaN
24+ vx_obs = NaN
25+ vy_obs = NaN
26+ vz_obs = NaN
27+ vel_obs = NaN
28+ thickness_obs = NaN
29+ surface_obs = NaN
30+
31+ end
32+ methods
33+ function self = adm1qn3inversion(varargin) % {{{
34+ switch nargin
35+ case 0
36+ self=setdefaultparameters(self);
37+ case 1
38+ self=structtoobj(adm1qn3inversion(),varargin{1});
39+ otherwise
40+ error('constructor not supported');
41+ end
42+ end % }}}
43+ function self = extrude(self,md) % {{{
44+ self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node');
45+ self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node');
46+ self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node');
47+ self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node');
48+ if numel(self.cost_functions_coefficients)>1,self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node');end;
49+ if numel(self.min_parameters)>1,self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node');end;
50+ if numel(self.max_parameters)>1,self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node');end;
51+ end % }}}
52+ function self = setdefaultparameters(self) % {{{
53+
54+
55+ %parameter to be inferred by control methods (only
56+ %drag and B are supported yet)
57+ self.control_parameters={'FrictionCoefficient'};
58+
59+ %Scaling factor for each control
60+ self.control_scaling_factors=1;
61+
62+ %number of iterations
63+ self.maxsteps=20;
64+ self.maxiter=40;
65+
66+ %several responses can be used:
67+ self.cost_functions=101;
68+
69+ %m1qn3 parameters
70+ self.dxmin = 0.1;
71+ self.gttol = 1e-4;
72+
73+ end % }}}
74+ function md = checkconsistency(self,md,solution,analyses) % {{{
75+
76+ %Early return
77+ if ~self.iscontrol, return; end
78+
79+ if ~IssmConfig('_HAVE_M1QN3_'),
80+ md = checkmessage(md,['M1QN3 has not been installed, ISSM needs to be reconfigured and recompiled with M1QN3']);
81+ end
82+ num_controls=numel(md.inversion.control_parameters);
83+ num_costfunc=size(md.inversion.cost_functions,2);
84+
85+ md = checkfield(md,'fieldname','inversion.iscontrol','values',[0 1]);
86+ md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',supportedcontrols());
87+ md = checkfield(md,'fieldname','inversion.control_scaling_factors','size',[1 num_controls],'>',0,'NaN',1,'Inf',1);
88+ md = checkfield(md,'fieldname','inversion.maxsteps','numel',1,'>=',0);
89+ md = checkfield(md,'fieldname','inversion.maxiter','numel',1,'>=',0);
90+ md = checkfield(md,'fieldname','inversion.dxmin','numel',1,'>',0);
91+ md = checkfield(md,'fieldname','inversion.gttol','numel',1,'>',0);
92+ md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',supportedcostfunctions());
93+ md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
94+ md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
95+ md = checkfield(md,'fieldname','inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]);
96+
97+ if strcmp(solution,'BalancethicknessSolution')
98+ md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
99+ md = checkfield(md,'fieldname','inversion.surface_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
100+ elseif strcmp(solution,'BalancethicknessSoftSolution')
101+ md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
102+ else
103+ md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
104+ if ~strcmp(domaintype(md.mesh),'2Dvertical'),
105+ md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
106+ end
107+ end
108+ end % }}}
109+ function disp(self) % {{{
110+ disp(sprintf(' adm1qn3inversion parameters:'));
111+ fielddisplay(self,'iscontrol','is inversion activated?');
112+ fielddisplay(self,'control_parameters','ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}');
113+ fielddisplay(self,'control_scaling_factors','order of magnitude of each control (useful for multi-parameter optimization)');
114+ fielddisplay(self,'maxsteps','maximum number of iterations (gradient computation)');
115+ fielddisplay(self,'maxiter','maximum number of Function evaluation (forward run)');
116+ fielddisplay(self,'dxmin','convergence criterion: two points less than dxmin from eachother (sup-norm) are considered identical');
117+ fielddisplay(self,'gttol','convergence criterion: ||g(X)||/||g(X0)|| (g(X0): gradient at initial guess X0)');
118+ fielddisplay(self,'cost_functions','indicate the type of response for each optimization step');
119+ fielddisplay(self,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter');
120+ fielddisplay(self,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex');
121+ fielddisplay(self,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex');
122+ fielddisplay(self,'vx_obs','observed velocity x component [m/yr]');
123+ fielddisplay(self,'vy_obs','observed velocity y component [m/yr]');
124+ fielddisplay(self,'vel_obs','observed velocity magnitude [m/yr]');
125+ fielddisplay(self,'thickness_obs','observed thickness [m]');
126+ fielddisplay(self,'surface_obs','observed surface elevation [m]');
127+ disp('Available cost functions:');
128+ disp(' 101: SurfaceAbsVelMisfit');
129+ disp(' 102: SurfaceRelVelMisfit');
130+ disp(' 103: SurfaceLogVelMisfit');
131+ disp(' 104: SurfaceLogVxVyMisfit');
132+ disp(' 105: SurfaceAverageVelMisfit');
133+ disp(' 201: ThicknessAbsMisfit');
134+ disp(' 501: DragCoefficientAbsGradient');
135+ disp(' 502: RheologyBbarAbsGradient');
136+ disp(' 503: ThicknessAbsGradient');
137+ end % }}}
138+ function marshall(self,prefix,md,fid) % {{{
139+
140+ yts=md.constants.yts;
141+
142+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','iscontrol','format','Boolean');
143+ WriteData(fid,prefix,'name','md.inversion.type','data',4,'format','Integer');
144+ if ~self.iscontrol, return; end
145+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','control_scaling_factors','format','DoubleMat','mattype',3);
146+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxsteps','format','Integer');
147+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxiter','format','Integer');
148+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','dxmin','format','Double');
149+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gttol','format','Double');
150+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
151+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3);
152+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3);
153+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts);
154+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts);
155+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts);
156+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vel_obs','format','DoubleMat','mattype',1,'scale',1./yts);
157+ if(numel(self.thickness_obs)==md.mesh.numberofelements),
158+ mattype=2;
159+ else
160+ mattype=1;
161+ end
162+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',mattype);
163+ WriteData(fid,prefix,'object',self,'class','inversion','fieldname','surface_obs','format','DoubleMat','mattype',mattype);
164+
165+ %process control parameters
166+ num_control_parameters=numel(self.control_parameters);
167+ WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray');
168+ WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer');
169+
170+ %process cost functions
171+ num_cost_functions=size(self.cost_functions,2);
172+ data=marshallcostfunctions(self.cost_functions);
173+ WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray');
174+ WriteData(fid,prefix,'data',num_cost_functions,'name','md.inversion.num_cost_functions','format','Integer');
175+ end % }}}
176+ function savemodeljs(self,fid,modelname) % {{{
177+
178+ writejsdouble(fid,[modelname '.inversion.iscontrol'],self.iscontrol);
179+ writejscellstring(fid,[modelname '.inversion.control_parameters'],self.control_parameters);
180+ writejsdouble(fid,[modelname '.inversion.control_scaling_factors'],self.control_scaling_factors);
181+ writejsdouble(fid,[modelname '.inversion.maxsteps'],self.maxsteps);
182+ writejsdouble(fid,[modelname '.inversion.maxiter'],self.maxiter);
183+ writejsdouble(fid,[modelname '.inversion.dxmin'],self.dxmin);
184+ writejsdouble(fid,[modelname '.inversion.gttol'],self.gttol);
185+ writejs2Darray(fid,[modelname '.inversion.cost_functions'],self.cost_functions);
186+ writejs2Darray(fid,[modelname '.inversion.cost_functions_coefficients'],self.cost_functions_coefficients);
187+ writejs1Darray(fid,[modelname '.inversion.min_parameters'],self.min_parameters);
188+ writejs1Darray(fid,[modelname '.inversion.max_parameters'],self.max_parameters);
189+ writejs1Darray(fid,[modelname '.inversion.vx_obs'],self.vx_obs);
190+ writejs1Darray(fid,[modelname '.inversion.vy_obs'],self.vy_obs);
191+ writejs1Darray(fid,[modelname '.inversion.vz_obs'],self.vz_obs);
192+ writejs1Darray(fid,[modelname '.inversion.vel_obs'],self.vel_obs);
193+ writejs1Darray(fid,[modelname '.inversion.thickness_obs'],self.thickness_obs);
194+ writejs1Darray(fid,[modelname '.inversion.surface_obs'],self.surface_obs);
195+
196+ end % }}}
197+ end
198+end
199Index: ../trunk-jpl/src/c/cores/WrapperCorePointerFromSolutionEnum.cpp
200===================================================================
201--- ../trunk-jpl/src/c/cores/WrapperCorePointerFromSolutionEnum.cpp (revision 22442)
202+++ ../trunk-jpl/src/c/cores/WrapperCorePointerFromSolutionEnum.cpp (revision 22443)
203@@ -46,7 +46,7 @@
204 case 1: solutioncore=controltao_core; break;
205 case 2: solutioncore=controlm1qn3_core; break;
206 case 3: solutioncore=controlvalidation_core; break;
207- case 4: solutioncore=controlad_core; break;
208+ case 4: solutioncore=controladm1qn3_core; break;
209 default: _error_("control type not supported");
210 }
211 }
212Index: ../trunk-jpl/src/c/cores/cores.h
213===================================================================
214--- ../trunk-jpl/src/c/cores/cores.h (revision 22442)
215+++ ../trunk-jpl/src/c/cores/cores.h (revision 22443)
216@@ -29,6 +29,7 @@
217 void control_core(FemModel* femmodel);
218 void controltao_core(FemModel* femmodel);
219 void controlm1qn3_core(FemModel* femmodel);
220+void controladm1qn3_core(FemModel* femmodel);
221 void controlad_core(FemModel* femmodel);
222 void controlvalidation_core(FemModel* femmodel);
223 void masstransport_core(FemModel* femmodel);
224Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp
225===================================================================
226--- ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 22442)
227+++ ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 22443)
228@@ -101,7 +101,7 @@
229 iomodel->FetchData(&control_scaling_factors,NULL,NULL,"md.inversion.control_scaling_factors");
230 parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_controls));
231 break;
232- case 4:/*M1QN3 AD*/
233+ case 4:/*AD M1QN3*/
234 parameters->AddObject(iomodel->CopyConstantObject("md.inversion.dxmin",InversionDxminEnum));
235 parameters->AddObject(iomodel->CopyConstantObject("md.inversion.gttol",InversionGttolEnum));
236 parameters->AddObject(iomodel->CopyConstantObject("md.inversion.maxsteps",InversionMaxstepsEnum));
237Index: ../trunk-jpl/src/c/Makefile.am
238===================================================================
239--- ../trunk-jpl/src/c/Makefile.am (revision 22442)
240+++ ../trunk-jpl/src/c/Makefile.am (revision 22443)
241@@ -282,6 +282,7 @@
242 ./cores/controltao_core.cpp\
243 ./cores/controlad_core.cpp\
244 ./cores/controlm1qn3_core.cpp\
245+ ./cores/controladm1qn3_core.cpp\
246 ./cores/controlvalidation_core.cpp\
247 ./cores/adjointstressbalance_core.cpp\
248 ./cores/adjointbalancethickness_core.cpp\
249Index: ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp
250===================================================================
251--- ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 22442)
252+++ ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 22443)
253@@ -1,4 +1,4 @@
254-/*!\file: controlm1qn3_core.cpp
255+/*!\file: controladm1qn3_core.cpp
256 * \brief: core of the control solution
257 */
258
259@@ -314,7 +314,7 @@
260 xDelete<IssmPDouble>(totalgradient);
261
262 }/*}}}*/
263-void controlm1qn3_core(FemModel* femmodel){/*{{{*/
264+void controladm1qn3_core(FemModel* femmodel){/*{{{*/
265
266 /*Intermediaries*/
267 long omode;
268@@ -466,5 +466,5 @@
269 }/*}}}*/
270
271 #else
272-void controlm1qn3_core(FemModel* femmodel){_error_("M1QN3 not installed");}
273+void controladm1qn3_core(FemModel* femmodel){_error_("M1QN3 not installed");}
274 #endif //_HAVE_M1QN3_
Note: See TracBrowser for help on using the repository browser.