source: issm/oecreview/Archive/24684-25833/ISSM-25255-25256.diff

Last change on this file was 25834, checked in by Mathieu Morlighem, 5 years ago

CHG: added 24684-25833
File size: 13.3 KB
RevLine 
[25834]1Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
2===================================================================
3--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 25255)
4+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 25256)
5@@ -45,6 +45,7 @@
6 #ifdef _HAVE_DAKOTA_
7 void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type);
8 void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);
9+ void InputScaleFromDakota(IssmDouble* distributed_values, IssmDouble* partition, int npart, int nt, int name);
10 #endif
11 void InputUpdateFromIoModel(int index, IoModel* iomodel);
12 void InputUpdateFromVector(IssmDouble* vector, int name, int type);
13Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
14===================================================================
15--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 25255)
16+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 25256)
17@@ -4706,6 +4706,101 @@
18 #endif
19
20 #ifdef _HAVE_DAKOTA_
21+void Penta::InputScaleFromDakota(IssmDouble* distributed_values, IssmDouble* partition, int npart, int nt, int name){/*{{{*/
22+
23+ int interp;
24+ int type;
25+
26+ /*Branch according to whether we have a transient or not input: */
27+ type=this->inputs2->GetInputObjectEnum(name);
28+ if(type==PentaInput2Enum){
29+ /*Figure out if we are P0 or P1 interpolation: */
30+ PentaInput2* pentainput = this->inputs2->GetPentaInput(name);
31+ PentaInput2* pentainput2 = this->inputs2->GetPentaInput(DummyEnum);
32+ interp=pentainput->GetInterpolation();
33+
34+ if (interp==P0Enum){
35+ /*Update the value if this element belongs to the partition: */
36+ if(partition[this->Sid()]!=-1){
37+ int lid=this->lid; pentainput->Serve(1,&lid);
38+ /*scale P0 value for this element, corresponding to the partition:*/
39+ IssmDouble value = pentainput->element_values[0];
40+ value*=distributed_values[(int)partition[this->Sid()]];
41+ pentainput2->SetInput(P0Enum,this->lid,value);
42+ }
43+ }
44+ else if (interp==P1Enum){
45+ IssmDouble values[NUMVERTICES];
46+ int lidlist[NUMVERTICES];
47+ this->GetVerticesLidList(&lidlist[0]);
48+ pentainput->Serve(NUMVERTICES,&lidlist[0]);
49+ for (int i=0;i<NUMVERTICES;i++){
50+ values[i]=pentainput->element_values[i];
51+ if(partition[this->vertices[i]->Sid()]!=-1) values[i]*=distributed_values[(int)partition[this->vertices[i]->Sid()]];
52+ }
53+ pentainput2->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
54+ }
55+ else _error_("Penta::InputScaleFromDakota error message: input interpolation " << EnumToStringx(interp) << " not supported yet!");
56+ }
57+ else if(type==TransientInput2Enum){
58+
59+
60+ IssmDouble* steps=NULL;
61+ int nsteps;
62+ TransientInput2* transientinput = NULL;
63+ TransientInput2* transientinput2 = NULL;
64+
65+ /*retrieve transient input:*/
66+ transientinput= this->inputs2->GetTransientInput(name);
67+ transientinput2= this->inputs2->GetTransientInput(DummyEnum);
68+
69+ /*retrieve time steps: */
70+ transientinput->GetAllTimes(&steps,&nsteps);
71+
72+ /*double check:*/
73+ if (nsteps!=nt && nt!=1) _error_("Penta:InputScaleFromDakota error message: transient input " << EnumToStringx(name) <<
74+ " should have the same number of time steps as the number of time values distributed by Dakota: " << nt << "\n");
75+
76+ /*needed to update inputs:*/
77+ int lidlist[NUMVERTICES];
78+ this->GetVerticesLidList(&lidlist[0]);
79+
80+ /*go through the transient inputs, and update:*/
81+ for (int i=0;i<nsteps;i++){
82+ PentaInput2* pentainput=transientinput->GetPentaInput(i);
83+ PentaInput2* pentainput2=transientinput2->GetPentaInput(i);
84+ interp=pentainput->GetInterpolation();
85+
86+ if (interp==P0Enum){
87+ /*Update the value if this element belongs to the partition: */
88+ if(partition[this->Sid()]!=-1){
89+ int lid=this->lid; pentainput->Serve(1,&lid);
90+ /*scale P0 value for this element, corresponding to the partition:*/
91+ IssmDouble value = pentainput->element_values[0];
92+ if(nt==1) value*=distributed_values[(int)partition[this->Sid()]]; //we scale all the time steps with the same distributed_value
93+ else value*=distributed_values[(int)partition[this->Sid()]*nsteps+i]; //we scale all the time steps with distributed value for each step
94+
95+ pentainput2->SetInput(P0Enum,this->lid,value);
96+ }
97+ }
98+ else if (interp==P1Enum){
99+ IssmDouble values[NUMVERTICES];
100+ pentainput->Serve(NUMVERTICES,&lidlist[0]);
101+ for (int j=0;j<NUMVERTICES;j++){
102+ values[j]=pentainput->element_values[j];
103+ if(partition[this->vertices[i]->Sid()]!=-1){
104+ if(nt==1) values[j]*=distributed_values[(int)partition[this->vertices[j]->Sid()]];//we scale all the time steps with the same distributed_value
105+ else values[j]*=distributed_values[(int)partition[this->vertices[j]->Sid()]*nsteps+i];//we scale all the time steps with distributed value for each step
106+ }
107+ }
108+ pentainput2->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
109+ }
110+ else _error_("Penta::InputScaleFromDakota error message: input interpolation " << EnumToStringx(interp) << " not supported yet!");
111+ }
112+ }
113+ else _error_("Penta::InputScaleFromDakota error message: input type " << EnumToStringx(name) << " not supported yet!");
114+}
115+/*}}}*/
116 void Penta::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
117
118 /*Check that name is an element input*/
119Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
120===================================================================
121--- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 25255)
122+++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 25256)
123@@ -199,8 +199,10 @@
124 void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input2* vx_input,Input2* vy_input,Input2* vz_input);
125
126 #ifdef _HAVE_DAKOTA_
127- void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type);
128- void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);
129+ void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type);
130+ void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);
131+ void InputScaleFromDakota(IssmDouble* distributed_values, IssmDouble* partition, int npart, int nt, int name);
132+
133 #endif
134
135 #ifdef _HAVE_GIA_
136Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
137===================================================================
138--- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 25255)
139+++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 25256)
140@@ -181,8 +181,10 @@
141 #endif
142
143 #ifdef _HAVE_DAKOTA_
144- void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type){_error_("not implemented yet");};
145- void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("not implemented yet");};
146+ void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type){_error_("not implemented yet");};
147+ void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("not implemented yet");};
148+ void InputScaleFromDakota(IssmDouble* distributed_values, IssmDouble* partition, int npart, int nt, int name){_error_("not implemented yet!");}
149+;
150 #endif
151 /*}}}*/
152 };
153Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h
154===================================================================
155--- ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 25255)
156+++ ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 25256)
157@@ -187,8 +187,9 @@
158 #endif
159
160 #ifdef _HAVE_DAKOTA_
161- void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("not implemented yet");};
162- void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type){_error_("not implemented yet");};
163+ void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("not implemented yet");};
164+ void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type){_error_("not implemented yet");};
165+ void InputScaleFromDakota(IssmDouble* distributed_values, IssmDouble* partition, int npart, int nt, int name){_error_("not implemented yet!");}
166 #endif
167 /*}}}*/
168 };
169Index: ../trunk-jpl/src/c/classes/Elements/Element.h
170===================================================================
171--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 25255)
172+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 25256)
173@@ -279,6 +279,7 @@
174 #ifdef _HAVE_DAKOTA_
175 virtual void InputUpdateFromMatrixDakota(IssmDouble* matrix, int rows, int ncols, int name, int type)=0;
176 virtual void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type)=0;
177+ virtual void InputScaleFromDakota(IssmDouble* distributed_values, IssmDouble* partition, int npart, int nt, int name)=0;
178 #endif
179 virtual void InputUpdateFromIoModel(int index, IoModel* iomodel)=0;
180 virtual void InputUpdateFromVector(IssmDouble* vector, int name, int type)=0;
181Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
182===================================================================
183--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25255)
184+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25256)
185@@ -6020,6 +6020,99 @@
186 #endif
187
188 #ifdef _HAVE_DAKOTA_
189+void Tria::InputScaleFromDakota(IssmDouble* distributed_values, IssmDouble* partition, int npart, int nt, int name){/*{{{*/
190+
191+ int interp;
192+ int type;
193+
194+ /*Branch according to whether we have a transient or not input: */
195+ type=this->inputs2->GetInputObjectEnum(name);
196+ if(type==TriaInput2Enum){
197+ /*Figure out if we are P0 or P1 interpolation: */
198+ TriaInput2* triainput = this->inputs2->GetTriaInput(name);
199+ TriaInput2* triainput2 = this->inputs2->GetTriaInput(DummyEnum);
200+ this->InputServe(triainput);
201+ interp=triainput->GetInterpolation();
202+
203+ if (interp==P0Enum){
204+ /*Update the value if this element belongs to the partition: */
205+ if(partition[this->Sid()]!=-1){
206+ /*scale P0 value for this element, corresponding to the partition:*/
207+ IssmDouble value = triainput->element_values[0];
208+ value*=distributed_values[(int)partition[this->Sid()]];
209+ triainput2->SetInput(P0Enum,this->lid,value);
210+ }
211+ }
212+ else if (interp==P1Enum){
213+ IssmDouble values[NUMVERTICES];
214+ int lidlist[NUMVERTICES];
215+ this->GetVerticesLidList(&lidlist[0]);
216+ for (int i=0;i<NUMVERTICES;i++){
217+ values[i]=triainput->element_values[i];
218+ if(partition[this->vertices[i]->Sid()]!=-1) values[i]*=distributed_values[(int)partition[this->vertices[i]->Sid()]];
219+ }
220+ triainput2->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
221+ }
222+ else _error_("Tria::InputScaleFromDakota error message: input interpolation " << EnumToStringx(interp) << " not supported yet!");
223+ }
224+ else if(type==TransientInput2Enum){
225+
226+
227+ IssmDouble* steps=NULL;
228+ int nsteps;
229+ TransientInput2* transientinput = NULL;
230+ TransientInput2* transientinput2 = NULL;
231+
232+ /*retrieve transient input:*/
233+ transientinput= this->inputs2->GetTransientInput(name);
234+ transientinput2= this->inputs2->GetTransientInput(DummyEnum);
235+
236+ /*retrieve time steps: */
237+ transientinput->GetAllTimes(&steps,&nsteps);
238+
239+ /*double check:*/
240+ if (nsteps!=nt && nt!=1) _error_("Tria:InputScaleFromDakota error message: transient input " << EnumToStringx(name) <<
241+ " should have the same number of time steps as the number of time values distributed by Dakota: " << nt << "\n");
242+
243+ /*needed to update inputs:*/
244+ int lidlist[NUMVERTICES];
245+ this->GetVerticesLidList(&lidlist[0]);
246+
247+ /*go through the transient inputs, and update:*/
248+ for (int i=0;i<nsteps;i++){
249+ TriaInput2* triainput=transientinput->GetTriaInput(i);
250+ TriaInput2* triainput2=transientinput2->GetTriaInput(i);
251+ this->InputServe(triainput);
252+ interp=triainput->GetInterpolation();
253+
254+ if (interp==P0Enum){
255+ /*Update the value if this element belongs to the partition: */
256+ if(partition[this->Sid()]!=-1){
257+ /*scale P0 value for this element, corresponding to the partition:*/
258+ IssmDouble value = triainput->element_values[0];
259+ if(nt==1) value*=distributed_values[(int)partition[this->Sid()]]; //we scale all the time steps with the same distributed_value
260+ else value*=distributed_values[(int)partition[this->Sid()]*nsteps+i]; //we scale all the time steps with distributed value for each step
261+
262+ triainput2->SetInput(P0Enum,this->lid,value);
263+ }
264+ }
265+ else if (interp==P1Enum){
266+ IssmDouble values[NUMVERTICES];
267+ for (int j=0;j<NUMVERTICES;j++){
268+ values[j]=triainput->element_values[j];
269+ if(partition[this->vertices[i]->Sid()]!=-1){
270+ if(nt==1) values[j]*=distributed_values[(int)partition[this->vertices[j]->Sid()]];//we scale all the time steps with the same distributed_value
271+ else values[j]*=distributed_values[(int)partition[this->vertices[j]->Sid()]*nsteps+i];//we scale all the time steps with distributed value for each step
272+ }
273+ }
274+ triainput2->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
275+ }
276+ else _error_("Tria::InputScaleFromDakota error message: input interpolation " << EnumToStringx(interp) << " not supported yet!");
277+ }
278+ }
279+ else _error_("Tria::InputScaleFromDakota error message: input type " << EnumToStringx(name) << " not supported yet!");
280+}
281+/*}}}*/
282 void Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
283
284 /*Check that name is an element input*/
Note: See TracBrowser for help on using the repository browser.