source: issm/oecreview/Archive/17984-18295/ISSM-18127-18128.diff@ 18296

Last change on this file since 18296 was 18296, checked in by Mathieu Morlighem, 11 years ago

Added 17984-18295
File size: 54.6 KB
RevLine 
[18296]1Index: ../trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp
2===================================================================
3--- ../trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp (revision 18127)
4+++ ../trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp (revision 18128)
5@@ -19,205 +19,226 @@
6 #include "./types.h"
7 #include "./isnan.h"
8
9-void BrentSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,void*),void* optargs){
10+void BrentSearch(IssmDouble** pJ,OptPars optpars,IssmDouble* X0,IssmDouble (*f)(IssmDouble*,void*),IssmDouble (*g)(IssmDouble**,IssmDouble*,void*),void* usr){
11
12 /* This routine is optimizing a given function using Brent's method
13 * (Golden or parabolic procedure)*/
14
15 /*Intermediary*/
16+ int iter;
17 IssmDouble si,gold,intervalgold,oldintervalgold;
18- IssmDouble parab_num,parab_den;
19- IssmDouble distance,cm_jump;
20+ IssmDouble parab_num,parab_den,distance;
21 IssmDouble fxmax,fxmin,fxbest;
22 IssmDouble fx,fx1,fx2;
23- IssmDouble xmax,xmin,xbest;
24- IssmDouble x,x1,x2,xm;
25+ IssmDouble x,x1,x2,xm,xbest;
26 IssmDouble tol1,tol2,seps;
27 IssmDouble tolerance = 1.e-4;
28- int maxiter ,iter;
29- bool loop= true,goldenflag;
30
31 /*Recover parameters:*/
32- xmin=optpars->xmin;
33- xmax=optpars->xmax;
34- maxiter=optpars->maxiter;
35- cm_jump=optpars->cm_jump;
36+ int nsteps = optpars.nsteps;
37+ int nsize = optpars.nsize;
38+ IssmDouble xmin = optpars.xmin;
39+ IssmDouble xmax = optpars.xmax;
40+ int *maxiter = optpars.maxiter;
41+ IssmDouble *cm_jump = optpars.cm_jump;
42
43- /*initialize counter and get response at the boundaries*/
44- iter=0;
45- fxmin = (*f)(xmin,optargs);
46- if (xIsNan<IssmDouble>(fxmin)) _error_("Function evaluation returned NaN");
47- cout<<setprecision(5);
48- if(VerboseControl()) _printf0_("\n");
49- if(VerboseControl()) _printf0_(" Iteration x f(x) Tolerance Procedure\n");
50- if(VerboseControl()) _printf0_("\n");
51- if(VerboseControl()) _printf0_(" N/A "<<setw(12)<<xmin<<" "<<setw(12)<<fxmin<<" N/A boundary\n");
52- fxmax = (*f)(xmax,optargs);
53- if (xIsNan<IssmDouble>(fxmax)) _error_("Function evaluation returned NaN");
54- if(VerboseControl()) _printf0_(" N/A "<<setw(12)<<xmax<<" "<<setw(12)<<fxmax<<" N/A boundary\n");
55+ /*Initialize gradient and controls*/
56+ IssmDouble* G = NULL;
57+ IssmDouble* J = xNew<IssmDouble>(nsteps);
58+ IssmDouble* X = xNew<IssmDouble>(nsize);
59
60- /*test if jump option activated and xmin==0*/
61- if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && (fxmax/fxmin)<cm_jump){
62- *psearch_scalar=xmax;
63- *pJ=fxmax;
64- return;
65- }
66+ /*start iterations*/
67+ for(int n=0;n<nsteps;n++){
68
69- /*initialize optimization variables*/
70- seps=sqrt(DBL_EPSILON); //precision of a IssmDouble
71- distance=0.0; //new_x=old_x + distance
72- gold=0.5*(3.0-sqrt(5.0)); //gold = 1 - golden ratio
73- intervalgold=0.0; //distance used by Golden procedure
74+ /*Reset some variables*/
75+ iter = 0;
76+ xmin = 0.;
77+ xmax = 1.;
78+ bool loop = true;
79+ cout<<setprecision(5);
80
81- /*1: initialize the values of the 4 x needed (x1,x2,x,xbest)*/
82- x1=xmin+gold*(xmax-xmin);
83- x2=x1;
84- xbest=x1;
85- x=xbest;
86+ /*Get current Gradient at xmin=0*/
87+ if(VerboseControl()) _printf0_("\n" << " step " << n+1 << "/" << nsteps << "\n");
88+ fxmin = (*g)(&G,X0,usr); if(xIsNan<IssmDouble>(fxmin)) _error_("Function evaluation returned NaN");
89+ if(VerboseControl()) _printf0_("\n");
90+ if(VerboseControl()) _printf0_(" Iteration x f(x) Tolerance Procedure\n");
91+ if(VerboseControl()) _printf0_("\n");
92+ if(VerboseControl()) _printf0_(" N/A "<<setw(12)<<xmin<<" "<<setw(12)<<fxmin<<" N/A boundary\n");
93+
94+ /*Get f(xmax)*/
95+ for(int i=0;i<nsize;i++) X[i]=X0[i]+xmax*G[i];
96+ fxmax = (*f)(X,usr); if (xIsNan<IssmDouble>(fxmax)) _error_("Function evaluation returned NaN");
97+ if(VerboseControl()) _printf0_(" N/A "<<setw(12)<<xmax<<" "<<setw(12)<<fxmax<<" N/A boundary\n");
98
99- /*2: call the function to be evaluated*/
100- fxbest = (*f)(x,optargs);
101- if(xIsNan<IssmDouble>(fxbest)) _error_("Function evaluation returned NaN");
102- iter=iter+1;
103+ /*test if jump option activated and xmin==0*/
104+ if(!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && (fxmax/fxmin)<cm_jump[n]){
105+ for(int i=0;i<nsize;i++) X0[i]=X0[i]+xmax*G[i];
106+ xDelete<IssmDouble>(G);
107+ J[n]=fxmax;
108+ continue;
109+ }
110
111- /*3: update the other variables*/
112- fx1=fxbest;
113- fx2=fxbest;
114- /*xm is always in the middle of a and b*/
115- xm=0.5*(xmin+xmax);
116- /*update tolerances*/
117- tol1=seps*sqrt(pow(xbest,2))+tolerance/3.0;
118- tol2=2.0*tol1;
119+ /*initialize optimization variables*/
120+ seps=sqrt(DBL_EPSILON); //precision of a IssmDouble
121+ distance=0.0; //new_x=old_x + distance
122+ gold=0.5*(3.0-sqrt(5.0)); //gold = 1 - golden ratio
123+ intervalgold=0.0; //distance used by Golden procedure
124
125- /*4: print result*/
126- if(VerboseControl())
127- _printf0_(" "<<setw(5)<<iter<<" "<<setw(12)<<xbest<<" "<<setw(12)<<fxbest<<" "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<" initial\n");
128- if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
129- if(VerboseControl()) _printf0_(" optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump << "\n");
130- loop=false;
131- }
132+ /*1: initialize the values of the 4 x needed (x1,x2,x,xbest)*/
133+ x1=xmin+gold*(xmax-xmin);
134+ x2=x1;
135+ xbest=x1;
136+ x=xbest;
137
138- while(loop){
139+ /*2: call the function to be evaluated*/
140+ for(int i=0;i<nsize;i++) X[i]=X0[i]+x*G[i];
141+ fxbest = (*f)(X,usr); if(xIsNan<IssmDouble>(fxbest)) _error_("Function evaluation returned NaN");
142+ iter++;
143
144- goldenflag=true;
145+ /*3: update the other variables*/
146+ fx1=fxbest;
147+ fx2=fxbest;
148+ /*xm is always in the middle of a and b*/
149+ xm=0.5*(xmin+xmax);
150+ /*update tolerances*/
151+ tol1=seps*sqrt(pow(xbest,2))+tolerance/3.0;
152+ tol2=2.0*tol1;
153
154- // Is a parabolic fit possible ?
155- if (sqrt(pow(intervalgold,2))>tol1){
156+ /*4: print result*/
157+ if(VerboseControl())
158+ _printf0_(" "<<setw(5)<<iter<<" "<<setw(12)<<xbest<<" "<<setw(12)<<fxbest<<" "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<" initial\n");
159+ if (!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && ((fxbest/fxmin)<cm_jump[n])){
160+ if(VerboseControl()) _printf0_(" optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump[n] << "\n");
161+ loop=false;
162+ }
163
164- // Yes, so fit parabola
165- goldenflag=false;
166- parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);;
167- parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1);
168+ while(loop){
169
170- //reverse p if necessary
171- if(parab_den>0.0){
172- parab_num=-parab_num;
173- }
174- parab_den=sqrt(pow(parab_den,2));
175- oldintervalgold=intervalgold;
176- intervalgold=distance;
177+ bool goldenflag=true;
178
179- // Is the parabola acceptable (we use seps here because in some configuration parab_num==parab_den*(xmax-xbest)
180- // and the result is not repeatable anymore
181- if (( sqrt(pow(parab_num,2)) < sqrt(pow(0.5*parab_den*oldintervalgold,2))) &&
182- (parab_num>parab_den*(xmin-xbest)+seps) &&
183- (parab_num<parab_den*(xmax-xbest)-seps)){
184+ // Is a parabolic fit possible ?
185+ if (sqrt(pow(intervalgold,2))>tol1){
186
187- // Yes, parabolic interpolation step
188- distance=parab_num/parab_den;
189- x=xbest+distance;
190+ // Yes, so fit parabola
191+ goldenflag=false;
192+ parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);;
193+ parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1);
194
195- // f must not be evaluated too close to min_x or max_x
196- if (((x-xmin)<tol2) || ((xmax-x)<tol2)){
197- if ((xm-xbest)<0.0) si=-1;
198- else si=1;
199- //compute new distance
200- distance=tol1*si;
201+ //reverse p if necessary
202+ if(parab_den>0.0){
203+ parab_num=-parab_num;
204 }
205+ parab_den=sqrt(pow(parab_den,2));
206+ oldintervalgold=intervalgold;
207+ intervalgold=distance;
208+
209+ // Is the parabola acceptable (we use seps here because in some configuration parab_num==parab_den*(xmax-xbest)
210+ // and the result is not repeatable anymore
211+ if (( sqrt(pow(parab_num,2)) < sqrt(pow(0.5*parab_den*oldintervalgold,2))) &&
212+ (parab_num>parab_den*(xmin-xbest)+seps) &&
213+ (parab_num<parab_den*(xmax-xbest)-seps)){
214+
215+ // Yes, parabolic interpolation step
216+ distance=parab_num/parab_den;
217+ x=xbest+distance;
218+
219+ // f must not be evaluated too close to min_x or max_x
220+ if (((x-xmin)<tol2) || ((xmax-x)<tol2)){
221+ if ((xm-xbest)<0.0) si=-1;
222+ else si=1;
223+ //compute new distance
224+ distance=tol1*si;
225+ }
226+ }
227+ else{
228+ // Not acceptable, must do a golden section step
229+ goldenflag=true;
230+ }
231 }
232- else{
233- // Not acceptable, must do a golden section step
234- goldenflag=true;
235- }
236- }
237
238- //Golden procedure
239- if(goldenflag){
240- // compute the new distance d
241- if(xbest>=xm){
242- intervalgold=xmin-xbest;
243+ //Golden procedure
244+ if(goldenflag){
245+ // compute the new distance d
246+ if(xbest>=xm){
247+ intervalgold=xmin-xbest;
248+ }
249+ else{
250+ intervalgold=xmax-xbest;
251+ }
252+ distance=gold*intervalgold;
253 }
254- else{
255- intervalgold=xmax-xbest;
256- }
257- distance=gold*intervalgold;
258- }
259
260- // The function must not be evaluated too close to xbest
261- if(distance<0) si=-1;
262- else si=1;
263- if(sqrt(pow(distance,2))>tol1) x=xbest+si*sqrt(pow(distance,2));
264- else x=xbest+si*tol1;
265+ // The function must not be evaluated too close to xbest
266+ if(distance<0) si=-1;
267+ else si=1;
268+ if(sqrt(pow(distance,2))>tol1) x=xbest+si*sqrt(pow(distance,2));
269+ else x=xbest+si*tol1;
270
271- //evaluate function on x
272- fx = (*f)(x,optargs);
273- if(xIsNan<IssmDouble>(fx)) _error_("Function evaluation returned NaN");
274- iter=iter+1;
275+ //evaluate function on x
276+ for(int i=0;i<nsize;i++) X[i]=X0[i]+x*G[i];
277+ fx = (*f)(X,usr); if(xIsNan<IssmDouble>(fx)) _error_("Function evaluation returned NaN");
278+ iter++;
279
280- // Update a, b, xm, x1, x2, tol1, tol2
281- if (fx<=fxbest){
282- if (x>=xbest) xmin=xbest;
283- else xmax=xbest;
284- x1=x2; fx1=fx2;
285- x2=xbest; fx2=fxbest;
286- xbest=x; fxbest=fx;
287- }
288- else{ // fx > fxbest
289- if (x<xbest) xmin=x;
290- else xmax=x;
291- if ((fx<=fx2) || (x2==xbest)){
292- x1=x2; fx1=fx2;
293- x2=x; fx2=fx;
294+ // Update a, b, xm, x1, x2, tol1, tol2
295+ if (fx<=fxbest){
296+ if (x>=xbest) xmin=xbest;
297+ else xmax=xbest;
298+ x1=x2; fx1=fx2;
299+ x2=xbest; fx2=fxbest;
300+ xbest=x; fxbest=fx;
301 }
302- else if ( (fx <= fx1) || (x1 == xbest) || (x1 == x2) ){
303- x1=x; fx1=fx;
304+ else{ // fx > fxbest
305+ if (x<xbest) xmin=x;
306+ else xmax=x;
307+ if ((fx<=fx2) || (x2==xbest)){
308+ x1=x2; fx1=fx2;
309+ x2=x; fx2=fx;
310+ }
311+ else if ( (fx <= fx1) || (x1 == xbest) || (x1 == x2) ){
312+ x1=x; fx1=fx;
313+ }
314 }
315- }
316- xm = 0.5*(xmin+xmax);
317- tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0;
318- tol2=2.0*tol1;
319- if(VerboseControl())
320- _printf0_(" "<<setw(5)<<iter<<" "<<setw(12)<<x<<" "<<setw(12)<<fx<<" "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<
321- " "<<(goldenflag?"golden":"parabolic")<<"\n");
322+ xm = 0.5*(xmin+xmax);
323+ tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0;
324+ tol2=2.0*tol1;
325+ if(VerboseControl())
326+ _printf0_(" "<<setw(5)<<iter<<" "<<setw(12)<<x<<" "<<setw(12)<<fx<<" "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<
327+ " "<<(goldenflag?"golden":"parabolic")<<"\n");
328
329- /*Stop the optimization?*/
330- if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(xmax-xmin))){
331- if(VerboseControl()) _printf0_(" optimization terminated: current x satisfies criteria 'tolx'=" << tolerance << "\n");
332- loop=false;
333+ /*Stop the optimization?*/
334+ if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(xmax-xmin))){
335+ if(VerboseControl()) _printf0_(" optimization terminated: current x satisfies criteria 'tolx'=" << tolerance << "\n");
336+ loop=false;
337+ }
338+ else if (iter>=maxiter[n]){
339+ if(VerboseControl()) _printf0_(" exiting: Maximum number of iterations has been exceeded ('maxiter'=" << maxiter[n] << ")\n");
340+ loop=false;
341+ }
342+ else if (!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && ((fxbest/fxmin)<cm_jump[n])){
343+ if(VerboseControl()) _printf0_(" optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump[n] << "\n");
344+ loop=false;
345+ }
346+ else{
347+ //continue
348+ loop=true;
349+ }
350+ }//end while
351+
352+ //Now, check that the value on the boundaries are not better than current fxbest
353+ if (fxbest>fxmin){
354+ xbest=optpars.xmin; fxbest=fxmin;
355 }
356- else if (iter>=maxiter){
357- if(VerboseControl()) _printf0_(" exiting: Maximum number of iterations has been exceeded ('maxiter'=" << maxiter << ")\n");
358- loop=false;
359+ if (fxbest>fxmax){
360+ xbest=optpars.xmax; fxbest=fxmax;
361 }
362- else if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
363- if(VerboseControl()) _printf0_(" optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump << "\n");
364- loop=false;
365- }
366- else{
367- //continue
368- loop=true;
369- }
370- }//end while
371
372- //Now, check that the value on the boundaries are not better than current fxbest
373- if (fxbest>fxmin){
374- xbest=optpars->xmin; fxbest=fxmin;
375+ /*Assign output pointers: */
376+ for(int i=0;i<nsize;i++) X0[i]=X0[i]+xbest*G[i];
377+ xDelete<IssmDouble>(G);
378+ J[n]=fxbest;
379 }
380- if (fxbest>fxmax){
381- xbest=optpars->xmax; fxbest=fxmax;
382- }
383-
384- /*Assign output pointers: */
385- *psearch_scalar=xbest;
386- *pJ=fxbest;
387+
388+ /*return*/
389+ xDelete<IssmDouble>(X);
390+ *pJ=J;
391 }
392Index: ../trunk-jpl/src/c/shared/Numerics/numerics.h
393===================================================================
394--- ../trunk-jpl/src/c/shared/Numerics/numerics.h (revision 18127)
395+++ ../trunk-jpl/src/c/shared/Numerics/numerics.h (revision 18128)
396@@ -29,7 +29,7 @@
397
398 int min(int a,int b);
399 int max(int a,int b);
400-void BrentSearch(IssmDouble *psearch_scalar,IssmDouble*pJ,OptPars*optpars,IssmDouble (*f)(IssmDouble,void*),void* optargs);
401+void BrentSearch(IssmDouble** pJ,OptPars optpars,IssmDouble* X0,IssmDouble (*f)(IssmDouble*,void*),IssmDouble (*g)(IssmDouble**,IssmDouble*,void*),void* usr);
402 void cross(IssmDouble *result,IssmDouble*vector1,IssmDouble*vector2);
403 void XZvectorsToCoordinateSystem(IssmDouble *T,IssmDouble*xzvectors);
404 int cubic(IssmDouble a, IssmDouble b, IssmDouble c, IssmDouble d, double X[3], int *num);
405Index: ../trunk-jpl/src/c/shared/Numerics/OptPars.h
406===================================================================
407--- ../trunk-jpl/src/c/shared/Numerics/OptPars.h (revision 18127)
408+++ ../trunk-jpl/src/c/shared/Numerics/OptPars.h (revision 18128)
409@@ -9,10 +9,12 @@
410
411 struct OptPars{
412
413- IssmDouble xmin;
414- IssmDouble xmax;
415- IssmDouble cm_jump;
416- int maxiter;
417+ IssmDouble xmin;
418+ IssmDouble xmax;
419+ IssmDouble *cm_jump;
420+ int* maxiter;
421+ int nsteps;
422+ int nsize;
423
424 };
425
426Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp
427===================================================================
428--- ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 18127)
429+++ ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 18128)
430@@ -15,10 +15,10 @@
431 int num_control_type;
432 int num_cm_responses;
433 int *control_type = NULL;
434+ int *maxiter = NULL;
435 int *cm_responses = NULL;
436 IssmDouble *cm_jump = NULL;
437 IssmDouble *optscal = NULL;
438- IssmDouble *maxiter = NULL;
439
440 /*retrieve some parameters: */
441 iomodel->Constant(&control_analysis,InversionIscontrolEnum);
442@@ -55,7 +55,7 @@
443 iomodel->FetchData(&maxiter,NULL,NULL,InversionMaxiterPerStepEnum);
444 parameters->AddObject(new DoubleMatParam(InversionGradientScalingEnum,optscal,nsteps,num_control_type));
445 parameters->AddObject(new DoubleVecParam(InversionStepThresholdEnum,cm_jump,nsteps));
446- parameters->AddObject(new DoubleVecParam(InversionMaxiterPerStepEnum,maxiter,nsteps));
447+ parameters->AddObject(new IntVecParam(InversionMaxiterPerStepEnum,maxiter,nsteps));
448 break;
449 case 1:/*TAO*/
450 parameters->AddObject(iomodel->CopyConstantObject(InversionFatolEnum));
451@@ -82,8 +82,8 @@
452
453 xDelete<int>(control_type);
454 xDelete<int>(cm_responses);
455+ xDelete<int>(maxiter);
456 iomodel->DeleteData(cm_jump,InversionStepThresholdEnum);
457 iomodel->DeleteData(optscal,InversionGradientScalingEnum);
458- iomodel->DeleteData(maxiter,InversionMaxiterPerStepEnum);
459 }
460 }
461Index: ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp
462===================================================================
463--- ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp (revision 18127)
464+++ ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp (revision 18128)
465@@ -6,18 +6,18 @@
466 #include "../../shared/shared.h"
467 #include "../../toolkits/toolkits.h"
468
469-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,IssmDouble* vector){
470+void SetControlInputsFromVectorx(FemModel* femmodel,IssmDouble* vector){
471
472 int num_controls;
473 int *control_type = NULL;
474
475 /*Retrieve some parameters*/
476- parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
477- parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
478+ femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
479+ femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
480
481 for(int i=0;i<num_controls;i++){
482- for(int j=0;j<elements->Size();j++){
483- Element* element=(Element*)elements->GetObjectByOffset(j);
484+ for(int j=0;j<femmodel->elements->Size();j++){
485+ Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
486 element->SetControlInputsFromVector(vector,control_type[i],i);
487 }
488 }
489@@ -25,14 +25,9 @@
490 xDelete<int>(control_type);
491 }
492
493-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,Vector<IssmDouble>* vector){
494+void SetControlInputsFromVectorx(FemModel* femmodel,Vector<IssmDouble>* vector){
495
496- IssmDouble* serial_vector=NULL;
497-
498- serial_vector=vector->ToMPISerial();
499-
500- SetControlInputsFromVectorx(elements,nodes, vertices, loads, materials, parameters,serial_vector);
501-
502- /*Free ressources:*/
503+ IssmDouble* serial_vector=vector->ToMPISerial();
504+ SetControlInputsFromVectorx(femmodel,serial_vector);
505 xDelete<IssmDouble>(serial_vector);
506 }
507Index: ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h
508===================================================================
509--- ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h (revision 18127)
510+++ ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h (revision 18128)
511@@ -7,7 +7,7 @@
512 #include "../../classes/classes.h"
513
514 /* local prototypes: */
515-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,Vector<IssmDouble>* vector);
516-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,IssmDouble* vector);
517+void SetControlInputsFromVectorx(FemModel* femmodel,Vector<IssmDouble>* vector);
518+void SetControlInputsFromVectorx(FemModel* femmodel,IssmDouble* vector);
519
520 #endif
521Index: ../trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp
522===================================================================
523--- ../trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 18127)
524+++ ../trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 18128)
525@@ -49,6 +49,8 @@
526 gradient->AXPY(gradient_list[i],1.0);
527 delete gradient_list[i];
528 }
529+ //gradient->Echo();
530+ //_error_("S");
531
532 /*Check that gradient is clean*/
533 norm_inf=gradient->Norm(NORM_INF);
534Index: ../trunk-jpl/src/c/cores/controlvalidation_core.cpp
535===================================================================
536--- ../trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 18127)
537+++ ../trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 18128)
538@@ -72,7 +72,7 @@
539 for(int i=0;i<n;i++) X[i] = X0[i] + alpha;
540
541 /*Calculate j(k+alpha delta k) */
542- SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
543+ SetControlInputsFromVectorx(femmodel,X);
544 solutioncore(femmodel);
545 femmodel->CostFunctionx(&j,NULL,NULL);
546
547Index: ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp
548===================================================================
549--- ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp (revision 18127)
550+++ ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp (revision 18128)
551@@ -105,7 +105,7 @@
552 }
553
554 /*Get solution*/
555- SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
556+ SetControlInputsFromVectorx(femmodel,X);
557 ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
558 femmodel->OutputControlsx(&femmodel->results);
559 femmodel->results->AddObject(new GenericExternalResult<double>(JEnum,f,1,0));
560@@ -130,10 +130,9 @@
561 int solution_type;
562 FemModel *femmodel = (FemModel*)dzs;
563
564- /*Recover responses*/
565- int num_responses;
566- int *responses = NULL;
567- femmodel->parameters->FindParam(&responses,&num_responses,InversionCostFunctionsEnum);
568+ /*Recover number of cost functions responses*/
569+ int num_responses;
570+ femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
571
572 /*Constrain input vector*/
573 IssmDouble *XL = NULL;
574@@ -146,7 +145,7 @@
575 }
576
577 /*Update control input*/
578- SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
579+ SetControlInputsFromVectorx(femmodel,X);
580
581 /*Recover some parameters*/
582 femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
583@@ -169,7 +168,6 @@
584
585 if(indic==0){
586 /*dry run, no gradient required*/
587- xDelete<int>(responses);
588 xDelete<IssmDouble>(XU);
589 xDelete<IssmDouble>(XL);
590 return;
591@@ -192,7 +190,6 @@
592 }
593
594 /*Clean-up and return*/
595- xDelete<int>(responses);
596 xDelete<IssmDouble>(XU);
597 xDelete<IssmDouble>(XL);
598 }
599Index: ../trunk-jpl/src/c/cores/control_core.cpp
600===================================================================
601--- ../trunk-jpl/src/c/cores/control_core.cpp (revision 18127)
602+++ ../trunk-jpl/src/c/cores/control_core.cpp (revision 18128)
603@@ -10,15 +10,19 @@
604 #include "../solutionsequences/solutionsequences.h"
605
606 /*Local prototypes*/
607-bool controlconvergence(IssmDouble J, IssmDouble tol_cm);
608-IssmDouble objectivefunction(IssmDouble search_scalar,void* optargs);
609+IssmDouble FormFunction(IssmDouble* X,void* usr);
610+IssmDouble FormFunctionGradient(IssmDouble** pG,IssmDouble* X,void* usr);
611+typedef struct {
612+ FemModel* femmodel;
613+ int nsize;
614+} AppCtx;
615
616 void control_core(FemModel* femmodel){
617
618 int i;
619
620 /*parameters: */
621- int num_controls;
622+ int num_controls,nsize;
623 int nsteps;
624 IssmDouble tol_cm;
625 int solution_type;
626@@ -26,11 +30,10 @@
627 bool dakota_analysis = false;
628
629 int *control_type = NULL;
630- IssmDouble *maxiter = NULL;
631+ int* maxiter = NULL;
632 IssmDouble *cm_jump = NULL;
633
634 /*intermediary: */
635- IssmDouble search_scalar = 1.;
636 OptPars optpars;
637
638 /*Solution and Adjoint core pointer*/
639@@ -60,37 +63,41 @@
640 if(VerboseControl()) _printf0_(" preparing initial solution\n");
641 if(isFS) solutioncore(femmodel);
642
643- /*Initialize cost function: */
644- J=xNew<IssmDouble>(nsteps);
645+ /*Get initial guess*/
646+ Vector<IssmDouble> *Xpetsc = NULL;
647+ GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
648+ IssmDouble* X0 = Xpetsc->ToMPISerial();
649+ Xpetsc->GetSize(&nsize);
650+ delete Xpetsc;
651
652 /*Initialize some of the BrentSearch arguments: */
653- optpars.xmin=0; optpars.xmax=1;
654+ optpars.xmin = 0;
655+ optpars.xmax = 1;
656+ optpars.nsteps = nsteps;
657+ optpars.nsize = nsize;
658+ optpars.maxiter = maxiter;
659+ optpars.cm_jump = cm_jump;
660
661- /*Start looping: */
662- for(int n=0;n<nsteps;n++){
663+ /*Initialize function argument*/
664+ AppCtx usr;
665+ usr.femmodel = femmodel;
666+ usr.nsize = nsize;
667
668- /*Display info*/
669- if(VerboseControl()) _printf0_("\n" << " control method step " << n+1 << "/" << nsteps << "\n");
670+ /*Call Brent optimization*/
671+ BrentSearch(&J,optpars,X0,&FormFunction,&FormFunctionGradient,(void*)&usr);
672
673-
674- /*In steady state inversion, compute new temperature field now*/
675- if(solution_type==SteadystateSolutionEnum) solutioncore(femmodel);
676-
677- if(VerboseControl()) _printf0_(" compute adjoint state:\n");
678- adjointcore(femmodel);
679- gradient_core(femmodel,n,search_scalar==0.);
680-
681- if(VerboseControl()) _printf0_(" optimizing along gradient direction\n");
682- optpars.maxiter=reCast<int,IssmDouble>(maxiter[n]); optpars.cm_jump=cm_jump[n];
683- BrentSearch(&search_scalar,J+n,&optpars,&objectivefunction,(void*)femmodel);
684-
685- if(VerboseControl()) _printf0_(" updating parameter using optimized search scalar\n"); //true means update save controls
686- InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,true);
687-
688- if(controlconvergence(J[n],tol_cm)) break;
689+ if(VerboseControl()) _printf0_(" preparing final solution\n");
690+ IssmDouble *XL = NULL;
691+ IssmDouble *XU = NULL;
692+ GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
693+ GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
694+ for(long i=0;i<nsize;i++){
695+ if(X0[i]>XU[i]) X0[i]=XU[i];
696+ if(X0[i]<XL[i]) X0[i]=XL[i];
697 }
698-
699- if(VerboseControl()) _printf0_(" preparing final solution\n");
700+ xDelete<IssmDouble>(XU);
701+ xDelete<IssmDouble>(XL);
702+ SetControlInputsFromVectorx(femmodel,X0);
703 femmodel->parameters->SetParam(true,SaveResultsEnum);
704 solutioncore(femmodel);
705
706@@ -110,84 +117,187 @@
707
708 /*Free ressources: */
709 xDelete<int>(control_type);
710- xDelete<IssmDouble>(maxiter);
711+ xDelete<int>(maxiter);
712 xDelete<IssmDouble>(cm_jump);
713 xDelete<IssmDouble>(J);
714+ xDelete<IssmDouble>(X0);
715 }
716-bool controlconvergence(IssmDouble J, IssmDouble tol_cm){
717
718- bool converged=false;
719+IssmDouble FormFunction(IssmDouble* X,void* usrvoid){
720
721- /*Has convergence been reached?*/
722- if (!xIsNan<IssmDouble>(tol_cm) && J<tol_cm){
723- converged=true;
724- if(VerboseConvergence()) _printf0_(" Convergence criterion reached: J = " << J << " < " << tol_cm);
725- }
726-
727- return converged;
728-}
729-
730-IssmDouble objectivefunction(IssmDouble search_scalar,void* optargs){
731-
732 /*output: */
733 IssmDouble J;
734
735 /*parameters: */
736- int solution_type,analysis_type;
737- bool isFS = false;
738+ int solution_type,analysis_type,num_responses;
739+ bool isFS = false;
740 bool conserve_loads = true;
741- FemModel *femmodel = (FemModel*)optargs;
742+ AppCtx* usr = (AppCtx*)usrvoid;
743+ FemModel *femmodel = usr->femmodel;
744+ int nsize = usr->nsize;
745
746 /*Recover parameters: */
747 femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
748 femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
749 femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
750+ femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
751
752- /*set analysis type to compute velocity: */
753+ /*Constrain input vector*/
754+ IssmDouble *XL = NULL;
755+ IssmDouble *XU = NULL;
756+ GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
757+ GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
758+ for(long i=0;i<nsize;i++){
759+ if(X[i]>XU[i]) X[i]=XU[i];
760+ if(X[i]<XL[i]) X[i]=XL[i];
761+ }
762+
763+ /*Update control input*/
764+ SetControlInputsFromVectorx(femmodel,X);
765+
766+ /*solve forward: */
767 switch(solution_type){
768 case SteadystateSolutionEnum:
769+ femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
770+ stressbalance_core(femmodel); //We need a 3D velocity!! (vz is required for the next thermal run)
771+ break;
772 case StressbalanceSolutionEnum:
773 femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
774+ solutionsequence_nonlinear(femmodel,conserve_loads);
775 break;
776 case BalancethicknessSolutionEnum:
777 femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
778+ solutionsequence_linear(femmodel);
779 break;
780 case BalancethicknessSoftSolutionEnum:
781- femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
782+ /*NOTHING*/
783 break;
784 case Balancethickness2SolutionEnum:
785 femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
786+ solutionsequence_linear(femmodel);
787 break;
788 default:
789 _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
790 }
791
792- /*update parameter according to scalar: */ //false means: do not save control
793- InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,false);
794+ /*Compute misfit for this velocity field.*/
795+ IssmDouble* Jlist = NULL;
796+ femmodel->CostFunctionx(&J,&Jlist,NULL);
797+ //_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<" | ");
798
799- /*Run stressbalance with updated inputs: */
800- if (solution_type==SteadystateSolutionEnum){
801- stressbalance_core(femmodel); //We need a 3D velocity!! (vz is required for the next thermal run)
802+ /*Retrieve objective functions independently*/
803+ //for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
804+ //_printf0_("\n");
805+
806+ /*Free ressources:*/
807+ xDelete<IssmDouble>(XU);
808+ xDelete<IssmDouble>(XL);
809+ xDelete<IssmDouble>(Jlist);
810+ return J;
811+}
812+IssmDouble FormFunctionGradient(IssmDouble** pG,IssmDouble* X,void* usrvoid){
813+
814+ /*output: */
815+ IssmDouble J;
816+
817+ /*parameters: */
818+ void (*adjointcore)(FemModel*)=NULL;
819+ int solution_type,analysis_type,num_responses,num_controls,numvertices;
820+ bool isFS = false;
821+ bool conserve_loads = true;
822+ IssmDouble *scalar_list = NULL;
823+ IssmDouble *Jlist = NULL;
824+ IssmDouble *G = NULL;
825+ IssmDouble *norm_list = NULL;
826+ AppCtx *usr = (AppCtx*)usrvoid;
827+ FemModel *femmodel = usr->femmodel;
828+ int nsize = usr->nsize;
829+
830+ /*Recover parameters: */
831+ femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
832+ femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
833+ femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
834+ femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
835+ femmodel->parameters->FindParam(&scalar_list,NULL,NULL,InversionGradientScalingEnum);
836+ femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); _assert_(num_controls);
837+ numvertices=femmodel->vertices->NumberOfVertices();
838+
839+ /*Constrain input vector*/
840+ IssmDouble *XL = NULL;
841+ IssmDouble *XU = NULL;
842+ GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
843+ GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
844+ for(long i=0;i<nsize;i++){
845+ if(X[i]>XU[i]) X[i]=XU[i];
846+ if(X[i]<XL[i]) X[i]=XL[i];
847 }
848- else if (solution_type==StressbalanceSolutionEnum){
849- solutionsequence_nonlinear(femmodel,conserve_loads);
850+
851+ /*Update control input*/
852+ SetControlInputsFromVectorx(femmodel,X);
853+
854+ /*Compute Adjoint*/
855+ AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
856+ adjointcore(femmodel);
857+
858+ /*Compute gradient*/
859+ Gradjx(&G,&norm_list,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
860+
861+ /*Compute scaling factor*/
862+ IssmDouble scalar = scalar_list[0]/norm_list[0];
863+ for(int i=1;i<num_controls;i++) scalar=min(scalar,scalar_list[i]/norm_list[i]);
864+
865+ /*Constrain Gradient*/
866+ for(int i=0;i<num_controls;i++){
867+ for(int j=0;j<numvertices;j++){
868+ G[i*numvertices+j] = scalar*G[i*numvertices+j];
869+ }
870 }
871- else if (solution_type==BalancethicknessSolutionEnum){
872- solutionsequence_linear(femmodel);
873+
874+ /*Needed for output results (FIXME: should be placed 6 lines below)*/
875+ ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
876+
877+ for(long i=0;i<nsize;i++){
878+ if(X[i]>=XU[i]) G[i]=0.;
879+ if(X[i]<=XL[i]) G[i]=0.;
880 }
881- else if (solution_type==Balancethickness2SolutionEnum){
882- solutionsequence_linear(femmodel);
883+
884+ /*solve forward: (FIXME: not needed actually...)*/
885+ switch(solution_type){
886+ case SteadystateSolutionEnum:
887+ femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
888+ stressbalance_core(femmodel); //We need a 3D velocity!! (vz is required for the next thermal run)
889+ break;
890+ case StressbalanceSolutionEnum:
891+ femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
892+ solutionsequence_nonlinear(femmodel,conserve_loads);
893+ break;
894+ case BalancethicknessSolutionEnum:
895+ femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
896+ solutionsequence_linear(femmodel);
897+ break;
898+ case BalancethicknessSoftSolutionEnum:
899+ /*NOTHING*/
900+ break;
901+ case Balancethickness2SolutionEnum:
902+ femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
903+ solutionsequence_linear(femmodel);
904+ break;
905+ default:
906+ _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
907 }
908- else if (solution_type==BalancethicknessSoftSolutionEnum){
909- /*Don't do anything*/
910- }
911- else{
912- _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
913- }
914
915 /*Compute misfit for this velocity field.*/
916- femmodel->CostFunctionx(&J,NULL,NULL);
917+ femmodel->CostFunctionx(&J,&Jlist,NULL);
918+ //_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<" | ");
919+ //for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
920+ //_printf0_("\n");
921
922- /*Free ressources:*/
923+ /*Clean-up and return*/
924+ xDelete<IssmDouble>(XU);
925+ xDelete<IssmDouble>(XL);
926+ xDelete<IssmDouble>(norm_list);
927+ xDelete<IssmDouble>(scalar_list);
928+ xDelete<IssmDouble>(Jlist);
929+ *pG = G;
930 return J;
931 }
932Index: ../trunk-jpl/src/c/cores/controltao_core.cpp
933===================================================================
934--- ../trunk-jpl/src/c/cores/controltao_core.cpp (revision 18127)
935+++ ../trunk-jpl/src/c/cores/controltao_core.cpp (revision 18128)
936@@ -92,7 +92,7 @@
937 TaoGetSolutionVector(tao,&X->pvector->vector);
938 G=new Vector<IssmDouble>(0); VecFree(&G->pvector->vector);
939 TaoGetGradientVector(tao,&G->pvector->vector);
940- SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
941+ SetControlInputsFromVectorx(femmodel,X);
942 ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
943 femmodel->OutputControlsx(&femmodel->results);
944 femmodel->results->AddObject(new GenericExternalResult<double*>(femmodel->results->Size()+1,JEnum,user.J,maxiter+3,1,1,0));
945@@ -112,11 +112,11 @@
946 TaoDestroy(&tao);
947 TaoFinalize();
948 }
949-int FormFunctionGradient(TaoSolver tao, Vec Xpetsc, IssmDouble *fcn,Vec G,void *userCtx){
950+int FormFunctionGradient(TaoSolver tao, Vec Xpetsc, IssmDouble *fcn,Vec G,void *uservoid){
951
952 /*Retreive arguments*/
953 int solution_type;
954- AppCtx *user = (AppCtx *)userCtx;
955+ AppCtx *user = (AppCtx *)uservoid;
956 FemModel *femmodel = user->femmodel;
957 Vector<IssmDouble> *gradient = NULL;
958 Vector<IssmDouble> *X = NULL;
959@@ -126,7 +126,7 @@
960
961 /*Set new variable*/
962 //VecView(X,PETSC_VIEWER_STDOUT_WORLD);
963- SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
964+ SetControlInputsFromVectorx(femmodel,X);
965 delete X;
966
967 /*Recover some parameters*/
968Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
969===================================================================
970--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18127)
971+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18128)
972@@ -983,7 +983,6 @@
973 this->inputs->AddInput(new IntInput(ApproximationEnum,reCast<int>(iomodel->Data(FlowequationElementEquationEnum)[index])));
974 }
975
976-
977 /*Control Inputs*/
978 if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
979 for(i=0;i<num_control_type;i++){
980Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
981===================================================================
982--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18127)
983+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18128)
984@@ -2850,6 +2850,7 @@
985 GradientIndexing(&vertexpidlist[0],control_index);
986
987 /*Get input (either in element or material)*/
988+ if(control_enum==MaterialsRheologyBbarEnum) control_enum=MaterialsRheologyBEnum;
989 Input* input=inputs->GetInput(control_enum);
990 if(!input) _error_("Input " << EnumToStringx(control_enum) << " not found in element");
991
992@@ -2864,10 +2865,21 @@
993 void Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
994
995 IssmDouble values[NUMVERTICES];
996- int vertexpidlist[NUMVERTICES];
997- Input *input = NULL;
998- Input *new_input = NULL;
999+ int vertexpidlist[NUMVERTICES],control_init;
1000+ Input *input = NULL;
1001+ Input *new_input = NULL;
1002
1003+ /*Specific case for depth averaged quantities*/
1004+ control_init=control_enum;
1005+ if(control_enum==MaterialsRheologyBbarEnum){
1006+ control_enum=MaterialsRheologyBEnum;
1007+ if(!IsOnBase()) return;
1008+ }
1009+ if(control_enum==DamageDbarEnum){
1010+ control_enum=DamageDEnum;
1011+ if(!IsOnBase()) return;
1012+ }
1013+
1014 /*Get out if this is not an element input*/
1015 if(!IsInput(control_enum)) return;
1016
1017@@ -2875,7 +2887,7 @@
1018 GradientIndexing(&vertexpidlist[0],control_index);
1019
1020 /*Get values on vertices*/
1021- for (int i=0;i<NUMVERTICES;i++){
1022+ for(int i=0;i<NUMVERTICES;i++){
1023 values[i]=vector[vertexpidlist[i]];
1024 }
1025 new_input = new PentaInput(control_enum,values,P1Enum);
1026@@ -2886,6 +2898,13 @@
1027 }
1028
1029 ((ControlInput*)input)->SetInput(new_input);
1030+
1031+ if(control_init==MaterialsRheologyBbarEnum){
1032+ this->InputExtrude(control_enum);
1033+ }
1034+ if(control_init==DamageDbarEnum){
1035+ this->InputExtrude(control_enum);
1036+ }
1037 }
1038 /*}}}*/
1039
1040Index: ../trunk-jpl/src/m/classes/inversion.m
1041===================================================================
1042--- ../trunk-jpl/src/m/classes/inversion.m (revision 18127)
1043+++ ../trunk-jpl/src/m/classes/inversion.m (revision 18128)
1044@@ -24,55 +24,55 @@
1045 thickness_obs = NaN
1046 end
1047 methods
1048- function createxml(obj,fid) % {{{
1049- fprintf(fid, '<!-- inversion -->\n');
1050-
1051- % inversion parameters
1052- fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="inversion parameters">','<section name="inversion" />');
1053- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="iscontrol" type="',class(obj.iscontrol),'" default="',convert2str(obj.iscontrol),'">',' <section name="inversion" />',' <help> is inversion activated? </help>',' </parameter>');
1054-
1055- % incompleteadjoing drop-down (0 or 1)
1056- fprintf(fid,'%s\n%s\n%s\n%s\n',' <parameter key ="incomplete_adjoint" type="alternative" optional="false">',' <section name="inversion" />',' <help> 1: linear viscosity, 0: non-linear viscosity </help>');
1057- fprintf(fid,'%s\n',' <option value="0" type="string" default="true"> </option>');
1058- fprintf(fid,'%s\n%s\n',' <option value="1" type="string" default="false"> </option>','</parameter>');
1059-
1060- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="control_parameters" type="',class(obj.control_parameters),'" default="',convert2str(obj.control_parameters),'">',' <section name="inversion" />',' <help> ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''} </help>',' </parameter>');
1061-
1062- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="nsteps" type="',class(obj.nsteps),'" default="',convert2str(obj.nsteps),'">',' <section name="inversion" />',' <help> number of optimization searches </help>',' </parameter>');
1063- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="cost_functions" type="',class(obj.cost_functions),'" default="',convert2str(obj.cost_functions),'">',' <section name="inversion" />',' <help> indicate the type of response for each optimization step </help>',' </parameter>');
1064- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="cost_functions_coefficients" type="',class(obj.cost_functions_coefficients),'" default="',convert2str(obj.cost_functions_coefficients),'">',' <section name="inversion" />',' <help> cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter </help>',' </parameter>');
1065-
1066- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="cost_function_threshold" type="',class(obj.cost_function_threshold),'" default="',convert2str(obj.cost_function_threshold),'">',' <section name="inversion" />',' <help> misfit convergence criterion. Default is 1%, NaN if not applied </help>',' </parameter>');
1067- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="maxiter_per_step" type="',class(obj.maxiter_per_step),'" default="',convert2str(obj.maxiter_per_step),'">',' <section name="inversion" />',' <help> maximum iterations during each optimization step </help>',' </parameter>');
1068- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="gradient_scaling" type="',class(obj.gradient_scaling),'" default="',convert2str(obj.gradient_scaling),'">',' <section name="inversion" />',' <help> scaling factor on gradient direction during optimization, for each optimization step </help>',' </parameter>');
1069-
1070- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="step_threshold" type="',class(obj.step_threshold),'" default="',convert2str(obj.step_threshold),'">',' <section name="inversion" />',' <help> decrease threshold for misfit, default is 30% </help>',' </parameter>');
1071- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="min_parameters" type="',class(obj.min_parameters),'" default="',convert2str(obj.min_parameters),'">',' <section name="inversion" />',' <help> absolute minimum acceptable value of the inversed parameter on each vertex </help>',' </parameter>');
1072- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="max_parameters" type="',class(obj.max_parameters),'" default="',convert2str(obj.max_parameters),'">',' <section name="inversion" />',' <help> absolute maximum acceptable value of the inversed parameter on each vertex </help>',' </parameter>');
1073-
1074- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="vx_obs" type="',class(obj.vx_obs),'" default="',convert2str(obj.vx_obs),'">',' <section name="inversion" />',' <help> observed velocity x component [m/yr] </help>',' </parameter>');
1075- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="vy_obs" type="',class(obj.vy_obs),'" default="',convert2str(obj.vy_obs),'">',' <section name="inversion" />',' <help> observed velocity y component [m/yr] </help>',' </parameter>');
1076- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="vel_obs" type="',class(obj.vel_obs),'" default="',convert2str(obj.vel_obs),'">',' <section name="inversion" />',' <help> observed velocity magnitude [m/yr] </help>',' </parameter>');
1077- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="thickness_obs" type="',class(obj.thickness_obs),'" default="',convert2str(obj.thickness_obs),'">',' <section name="inversion" />',' <help> observed thickness [m]) </help>',' </parameter>');
1078-
1079- fprintf(fid,'%s\n%s\n','</frame>');
1080-
1081- fprintf(fid,'%s\n%s\n%s\n','<frame key="2" label="Available cost functions">','<section name="inversion" />');
1082- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="SurfaceAbsVelMisfit" type="','string','" default="','101','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1083- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="SurfaceRelVelMisfit" type="','string','" default="','102','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1084- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="SurfaceLogVelMisfit" type="','string','" default="','103','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1085-
1086- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="SurfaceLogVxVyMisfit" type="','string','" default="','104','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1087- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="SurfaceAverageVelMisfit" type="','string','" default="','105','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1088- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="ThicknessAbsMisfit" type="','string','" default="','106','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1089-
1090- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="DragCoefficientAbsGradient" type="','string','" default="','107','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1091- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="RheologyBbarAbsGradient" type="','string','" default="','108','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1092- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="ThicknessAbsGradient" type="','string','" default="','109','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1093-
1094- fprintf(fid,'%s\n%s\n','</frame>');
1095-
1096- end % }}}
1097+ function createxml(obj,fid) % {{{
1098+ fprintf(fid, '<!-- inversion -->\n');
1099+
1100+ % inversion parameters
1101+ fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="inversion parameters">','<section name="inversion" />');
1102+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="iscontrol" type="',class(obj.iscontrol),'" default="',convert2str(obj.iscontrol),'">',' <section name="inversion" />',' <help> is inversion activated? </help>',' </parameter>');
1103+
1104+ % incompleteadjoing drop-down (0 or 1)
1105+ fprintf(fid,'%s\n%s\n%s\n%s\n',' <parameter key ="incomplete_adjoint" type="alternative" optional="false">',' <section name="inversion" />',' <help> 1: linear viscosity, 0: non-linear viscosity </help>');
1106+ fprintf(fid,'%s\n',' <option value="0" type="string" default="true"> </option>');
1107+ fprintf(fid,'%s\n%s\n',' <option value="1" type="string" default="false"> </option>','</parameter>');
1108+
1109+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="control_parameters" type="',class(obj.control_parameters),'" default="',convert2str(obj.control_parameters),'">',' <section name="inversion" />',' <help> ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''} </help>',' </parameter>');
1110+
1111+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="nsteps" type="',class(obj.nsteps),'" default="',convert2str(obj.nsteps),'">',' <section name="inversion" />',' <help> number of optimization searches </help>',' </parameter>');
1112+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="cost_functions" type="',class(obj.cost_functions),'" default="',convert2str(obj.cost_functions),'">',' <section name="inversion" />',' <help> indicate the type of response for each optimization step </help>',' </parameter>');
1113+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="cost_functions_coefficients" type="',class(obj.cost_functions_coefficients),'" default="',convert2str(obj.cost_functions_coefficients),'">',' <section name="inversion" />',' <help> cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter </help>',' </parameter>');
1114+
1115+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="cost_function_threshold" type="',class(obj.cost_function_threshold),'" default="',convert2str(obj.cost_function_threshold),'">',' <section name="inversion" />',' <help> misfit convergence criterion. Default is 1%, NaN if not applied </help>',' </parameter>');
1116+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="maxiter_per_step" type="',class(obj.maxiter_per_step),'" default="',convert2str(obj.maxiter_per_step),'">',' <section name="inversion" />',' <help> maximum iterations during each optimization step </help>',' </parameter>');
1117+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="gradient_scaling" type="',class(obj.gradient_scaling),'" default="',convert2str(obj.gradient_scaling),'">',' <section name="inversion" />',' <help> scaling factor on gradient direction during optimization, for each optimization step </help>',' </parameter>');
1118+
1119+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="step_threshold" type="',class(obj.step_threshold),'" default="',convert2str(obj.step_threshold),'">',' <section name="inversion" />',' <help> decrease threshold for misfit, default is 30% </help>',' </parameter>');
1120+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="min_parameters" type="',class(obj.min_parameters),'" default="',convert2str(obj.min_parameters),'">',' <section name="inversion" />',' <help> absolute minimum acceptable value of the inversed parameter on each vertex </help>',' </parameter>');
1121+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="max_parameters" type="',class(obj.max_parameters),'" default="',convert2str(obj.max_parameters),'">',' <section name="inversion" />',' <help> absolute maximum acceptable value of the inversed parameter on each vertex </help>',' </parameter>');
1122+
1123+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="vx_obs" type="',class(obj.vx_obs),'" default="',convert2str(obj.vx_obs),'">',' <section name="inversion" />',' <help> observed velocity x component [m/yr] </help>',' </parameter>');
1124+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="vy_obs" type="',class(obj.vy_obs),'" default="',convert2str(obj.vy_obs),'">',' <section name="inversion" />',' <help> observed velocity y component [m/yr] </help>',' </parameter>');
1125+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="vel_obs" type="',class(obj.vel_obs),'" default="',convert2str(obj.vel_obs),'">',' <section name="inversion" />',' <help> observed velocity magnitude [m/yr] </help>',' </parameter>');
1126+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="thickness_obs" type="',class(obj.thickness_obs),'" default="',convert2str(obj.thickness_obs),'">',' <section name="inversion" />',' <help> observed thickness [m]) </help>',' </parameter>');
1127+
1128+ fprintf(fid,'%s\n%s\n','</frame>');
1129+
1130+ fprintf(fid,'%s\n%s\n%s\n','<frame key="2" label="Available cost functions">','<section name="inversion" />');
1131+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="SurfaceAbsVelMisfit" type="','string','" default="','101','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1132+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="SurfaceRelVelMisfit" type="','string','" default="','102','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1133+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="SurfaceLogVelMisfit" type="','string','" default="','103','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1134+
1135+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="SurfaceLogVxVyMisfit" type="','string','" default="','104','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1136+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="SurfaceAverageVelMisfit" type="','string','" default="','105','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1137+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="ThicknessAbsMisfit" type="','string','" default="','106','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1138+
1139+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="DragCoefficientAbsGradient" type="','string','" default="','107','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1140+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="RheologyBbarAbsGradient" type="','string','" default="','108','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1141+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',' <parameter key ="ThicknessAbsGradient" type="','string','" default="','109','">',' <section name="inversion" />',' <help> </help>',' </parameter>');
1142+
1143+ fprintf(fid,'%s\n%s\n','</frame>');
1144+
1145+ end % }}}
1146 function obj = inversion(varargin) % {{{
1147 switch nargin
1148 case 0
1149@@ -195,7 +195,7 @@
1150 WriteData(fid,'object',obj,'fieldname','incomplete_adjoint','format','Boolean');
1151 if ~obj.iscontrol, return; end
1152 WriteData(fid,'object',obj,'fieldname','nsteps','format','Integer');
1153- WriteData(fid,'object',obj,'fieldname','maxiter_per_step','format','DoubleMat','mattype',3);
1154+ WriteData(fid,'object',obj,'fieldname','maxiter_per_step','format','IntMat','mattype',3);
1155 WriteData(fid,'object',obj,'fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
1156 WriteData(fid,'object',obj,'fieldname','gradient_scaling','format','DoubleMat','mattype',3);
1157 WriteData(fid,'object',obj,'fieldname','cost_function_threshold','format','Double');
Note: See TracBrowser for help on using the repository browser.