source: issm/oecreview/Archive/23185-23389/ISSM-23305-23306.diff@ 23390

Last change on this file since 23390 was 23390, checked in by Mathieu Morlighem, 6 years ago

CHG: added Archive/23185-23389
File size: 9.0 KB
RevLine 
[23390]1Index: ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp
2===================================================================
3--- ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp (revision 23305)
4+++ ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp (revision 23306)
5@@ -10,7 +10,7 @@
6 #include "../modules/modules.h"
7 #include "../solutionsequences/solutionsequences.h"
8
9-#if defined (_HAVE_M1QN3_) && !defined(_HAVE_AD_)
10+#if defined (_HAVE_M1QN3_)
11 /*m1qn3 prototypes {{{*/
12 extern "C" void *ctonbe_; // DIS mode : Conversion
13 extern "C" void *ctcabe_; // DIS mode : Conversion
14@@ -27,8 +27,8 @@
15
16 /*Use struct to provide arguments*/
17 typedef struct{
18- FemModel * femmodel;
19- IssmDouble* Jlist;
20+ FemModel * femmodel;
21+ IssmPDouble* Jlist;
22 int M;
23 int N;
24 int* i;
25@@ -43,8 +43,8 @@
26 int maxsteps,maxiter;
27 int intn,numberofvertices,num_controls,num_cost_functions,solution_type;
28 IssmDouble *scaling_factors = NULL;
29- IssmDouble *X = NULL;
30- IssmDouble *G = NULL;
31+ double *X = NULL;
32+ double *G = NULL;
33
34 /*Recover some parameters*/
35 femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
36@@ -52,8 +52,8 @@
37 femmodel->parameters->FindParam(&num_cost_functions,InversionNumCostFunctionsEnum);
38 femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum);
39 femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum);
40- femmodel->parameters->FindParam(&dxmin,InversionDxminEnum);
41- femmodel->parameters->FindParam(&gttol,InversionGttolEnum);
42+ femmodel->parameters->FindParamAndMakePassive(&dxmin,InversionDxminEnum);
43+ femmodel->parameters->FindParamAndMakePassive(&gttol,InversionGttolEnum);
44 femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
45 femmodel->parameters->SetParam(false,SaveResultsEnum);
46 numberofvertices=femmodel->vertices->NumberOfVertices();
47@@ -77,8 +77,8 @@
48 long nsim = long(maxiter);/*Maximum number of function calls*/
49
50 /*Get initial guess*/
51- Vector<IssmDouble> *Xpetsc = NULL;
52- GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
53+ Vector<double> *Xpetsc = NULL;
54+ GetPassiveVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
55 X = Xpetsc->ToMPISerial();
56 Xpetsc->GetSize(&intn);
57 delete Xpetsc;
58@@ -92,7 +92,7 @@
59 for(int c=0;c<num_controls;c++){
60 for(int i=0;i<numberofvertices;i++){
61 int index = numberofvertices*c+i;
62- X[index] = X[index]/scaling_factors[c];
63+ X[index] = X[index]/reCast<double>(scaling_factors[c]);
64 }
65 }
66
67@@ -111,7 +111,7 @@
68 mystruct.femmodel = femmodel;
69 mystruct.M = maxiter;
70 mystruct.N = num_cost_functions+1;
71- mystruct.Jlist = xNewZeroInit<IssmDouble>(mystruct.M*mystruct.N);
72+ mystruct.Jlist = xNewZeroInit<IssmPDouble>(mystruct.M*mystruct.N);
73 mystruct.i = xNewZeroInit<int>(1);
74
75 /*Initialize Gradient and cost function of M1QN3*/
76@@ -140,20 +140,36 @@
77 }
78
79 /*Constrain solution vector*/
80- IssmDouble *XL = NULL;
81- IssmDouble *XU = NULL;
82- GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
83- GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
84+ double *XL = NULL;
85+ double *XU = NULL;
86+ GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
87+ GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
88 for(int c=0;c<num_controls;c++){
89 for(int i=0;i<numberofvertices;i++){
90 int index = numberofvertices*c+i;
91- X[index] = X[index]*scaling_factors[c];
92+ X[index] = X[index]*reCast<double>(scaling_factors[c]);
93 if(X[index]>XU[index]) X[index]=XU[index];
94 if(X[index]<XL[index]) X[index]=XL[index];
95 }
96 }
97+
98+ /*Set X as our new control (need to recast)*/
99+ #ifdef _HAVE_AD_
100+ IssmDouble* aX=xNew<IssmDouble>(intn);
101+ IssmDouble* aG=xNew<IssmDouble>(intn);
102+ for(int i=0;i<intn;i++) {
103+ aX[i] = reCast<IssmDouble>(X[i]);
104+ aG[i] = reCast<IssmDouble>(G[i]);
105+ }
106+ ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
107+ SetControlInputsFromVectorx(femmodel,aX);
108+ xDelete(aX);
109+ xDelete(aG);
110+ #else
111 SetControlInputsFromVectorx(femmodel,X);
112 ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
113+ #endif
114+
115 femmodel->OutputControlsx(&femmodel->results);
116 femmodel->results->AddObject(new GenericExternalResult<double*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0));
117
118@@ -170,7 +186,7 @@
119 xDelete<double>(dz);
120 xDelete<double>(XU);
121 xDelete<double>(XL);
122- xDelete<double>(scaling_factors);
123+ xDelete<IssmDouble>(scaling_factors);
124 xDelete<double>(mystruct.Jlist);
125 xDelete<int>(mystruct.i);
126 }/*}}}*/
127@@ -179,7 +195,7 @@
128 /*Recover Arguments*/
129 m1qn3_struct *input_struct = (m1qn3_struct*)dzs;
130 FemModel *femmodel = input_struct->femmodel;
131- IssmDouble *Jlist = input_struct->Jlist;
132+ IssmPDouble *Jlist = input_struct->Jlist;
133 int JlistM = input_struct->M;
134 int JlistN = input_struct->N;
135 int *Jlisti = input_struct->i;
136@@ -194,19 +210,31 @@
137 numberofvertices=femmodel->vertices->NumberOfVertices();
138
139 /*Constrain input vector and update controls*/
140- IssmDouble *XL = NULL;
141- IssmDouble *XU = NULL;
142+ double *XL = NULL;
143+ double *XU = NULL;
144+ #ifdef _HAVE_AD_
145+ GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
146+ GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
147+ #else
148 GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
149 GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
150+ #endif
151 for(int c=0;c<num_controls;c++){
152 for(int i=0;i<numberofvertices;i++){
153 int index = numberofvertices*c+i;
154- X[index] = X[index]*scaling_factors[c];
155+ X[index] = X[index]*reCast<double>(scaling_factors[c]);
156 if(X[index]>XU[index]) X[index]=XU[index];
157 if(X[index]<XL[index]) X[index]=XL[index];
158 }
159 }
160+ #ifdef _HAVE_AD_
161+ IssmDouble* aX=xNew<IssmDouble>(*n);
162+ for(int i=0;i<*n;i++) aX[i] = reCast<IssmDouble>(X[i]);
163+ SetControlInputsFromVectorx(femmodel,aX);
164+ xDelete(aX);
165+ #else
166 SetControlInputsFromVectorx(femmodel,X);
167+ #endif
168
169 /*Compute solution and adjoint*/
170 void (*solutioncore)(FemModel*)=NULL;
171@@ -220,11 +248,13 @@
172
173 /*Compute objective function*/
174 IssmDouble* Jtemp = NULL;
175- femmodel->CostFunctionx(pf,&Jtemp,NULL);
176+ IssmDouble J;
177+ femmodel->CostFunctionx(&J,&Jtemp,NULL);
178+ *pf = reCast<double>(J);
179 _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<" | ");
180
181 /*Record cost function values and delete Jtemp*/
182- for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = Jtemp[i];
183+ for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = reCast<double>(Jtemp[i]);
184 Jlist[(*Jlisti)*JlistN+num_responses] = *pf;
185 xDelete<IssmDouble>(Jtemp);
186
187@@ -237,8 +267,8 @@
188 _printf0_("\n");
189
190 *Jlisti = (*Jlisti) +1;
191- xDelete<IssmDouble>(XU);
192- xDelete<IssmDouble>(XL);
193+ xDelete<double>(XU);
194+ xDelete<double>(XL);
195 return;
196 }
197
198@@ -249,7 +279,7 @@
199 /*Compute gradient*/
200 IssmDouble* G2 = NULL;
201 Gradjx(&G2,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
202- for(long i=0;i<*n;i++) G[i] = -G2[i];
203+ for(long i=0;i<*n;i++) G[i] = -reCast<double>(G2[i]);
204 xDelete<IssmDouble>(G2);
205
206 /*Constrain Gradient*/
207@@ -259,8 +289,8 @@
208 int index = numberofvertices*c+i;
209 if(X[index]>=XU[index]) G[index]=0.;
210 if(X[index]<=XL[index]) G[index]=0.;
211- G[index] = G[index]*scaling_factors[c];
212- X[index] = X[index]/scaling_factors[c];
213+ G[index] = G[index]*reCast<double>(scaling_factors[c]);
214+ X[index] = X[index]/reCast<double>(scaling_factors[c]);
215 Gnorm += G[index]*G[index];
216 }
217 }
218@@ -273,11 +303,11 @@
219
220 /*Clean-up and return*/
221 *Jlisti = (*Jlisti) +1;
222- xDelete<IssmDouble>(XU);
223- xDelete<IssmDouble>(XL);
224+ xDelete<double>(XU);
225+ xDelete<double>(XL);
226 xDelete<IssmDouble>(scaling_factors);
227 }/*}}}*/
228
229 #else
230-void controlm1qn3_core(FemModel* femmodel){_error_("M1QN3 not installed (or you might have turned AD on)");}
231+void controlm1qn3_core(FemModel* femmodel){_error_("M1QN3 not installed");}
232 #endif //_HAVE_M1QN3_
Note: See TracBrowser for help on using the repository browser.