source: issm/oecreview/Archive/24684-25833/ISSM-25488-25489.diff@ 25834

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

CHG: added 24684-25833
File size: 14.7 KB
RevLine 
[25834]1Index: ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 25488)
4+++ ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 25489)
5@@ -35,6 +35,44 @@
6 /*Assign output pointers:*/
7 *pvector=vector;
8 } /*}}}*/
9+void GetVectoronBaseFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type){ /*{{{*/
10+
11+ int i, domaintype;
12+ Vector<IssmDouble>* vector=NULL;
13+
14+ switch(type){
15+ case ElementSIdEnum:
16+ vector=new Vector<IssmDouble>(femmodel->elements->NumberOfElements());
17+ break;
18+ case VertexPIdEnum: case VertexSIdEnum:
19+ vector=new Vector<IssmDouble>(femmodel->vertices->NumberOfVertices());
20+ break;
21+ case NodesEnum:case NodeSIdEnum:
22+ vector=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes());
23+ break;
24+ default:
25+ _error_("vector type: " << EnumToStringx(type) << " not supported yet!");
26+ }
27+
28+ /*Look up in elements*/
29+ for(i=0;i<femmodel->elements->Size();i++){
30+ Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
31+ element->FindParam(&domaintype,DomainTypeEnum);
32+ switch(domaintype){
33+ case Domain2DhorizontalEnum:
34+ element->GetVectorFromInputs(vector,name,type);
35+ break;
36+ case Domain3DEnum:
37+ if(!element->IsOnBase()) continue;
38+ element->GetVectorFromInputs(vector,name,type);
39+ break;
40+ default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
41+ }
42+ }
43+ vector->Assemble();
44+ /*Assign output pointers:*/
45+ *pvector=vector;
46+} /*}}}*/
47 void GetVectorFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type,IssmDouble time){/*{{{*/
48
49 int i;
50@@ -78,18 +116,35 @@
51 /*Assign output pointers:*/
52 *pvector=vector;
53 }/*}}}*/
54+void GetVectoronBaseFromInputsx(IssmDouble** pvector,FemModel* femmodel,int name, int type){/*{{{*/
55+
56+ /*output: */
57+ IssmDouble* vector=NULL;
58+
59+ /*intermediary: */
60+ Vector<IssmDouble>* vec_vector=NULL;
61+
62+ GetVectoronBaseFromInputsx(&vec_vector,femmodel,name,type);
63+ vector=vec_vector->ToMPISerial();
64+
65+ /*Free ressources:*/
66+ delete vec_vector;
67+
68+ /*Assign output pointers:*/
69+ *pvector=vector;
70+}/*}}}*/
71 void GetVectorFromInputsx(IssmDouble** pvector,int* pvector_size, FemModel* femmodel,int name){ /*{{{*/
72
73 int interpolation_type;
74- /*this one is special: we don't specify the type, but let the nature of the inputs dictace.
75+ /*this one is special: we don't specify the type, but let the nature of the inputs dictace.
76 * P0 -> ElementSIdEnum, P1 ->VertexSIdEnum: */
77
78 /*We go find the input of the first element, and query its interpolation type: */
79 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(0));
80- Input* input=element->GetInput(name);
81+ Input* input=element->GetInput(name);
82 if (!input) _error_("could not find input: " << name);
83
84- interpolation_type=input->GetInputInterpolationType();
85+ interpolation_type=input->GetInputInterpolationType();
86 if(interpolation_type==P0Enum){
87 *pvector_size=femmodel->elements->NumberOfElements();
88 GetVectorFromInputsx(pvector,femmodel,name, ElementSIdEnum);
89Index: ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h
90===================================================================
91--- ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h (revision 25488)
92+++ ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h (revision 25489)
93@@ -1,5 +1,5 @@
94 /*!\file: GetVectorFromInputsx.h
95- */
96+ */
97
98 #ifndef _GETVECTORFROMINPUTSXX_H
99 #define _GETVECTORFROMINPUTSXX_H
100@@ -7,9 +7,12 @@
101 #include "../../classes/classes.h"
102
103 /* local prototypes: */
104-void GetVectorFromInputsx( Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type);
105-void GetVectorFromInputsx( IssmDouble** pvector,FemModel* femmodel,int name,int type);
106-void GetVectorFromInputsx(IssmDouble** pvector,int* pvector_size, FemModel* femmodel,int name);
107-void GetVectorFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type,IssmDouble time);
108+void GetVectorFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type);
109+void GetVectoronBaseFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type);
110+void GetVectorFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type,IssmDouble time);
111+void GetVectorFromInputsx(IssmDouble** pvector,FemModel* femmodel,int name,int type);
112+void GetVectoronBaseFromInputsx(IssmDouble** pvector,FemModel* femmodel,int name,int type);
113+void GetVectorFromInputsx(IssmDouble** pvector,int* pvector_size, FemModel* femmodel,int name);
114
115+
116 #endif /* _GETVECTORFROMINPUTSXX_H */
117Index: ../trunk-jpl/src/c/cores/hydrology_core.cpp
118===================================================================
119--- ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 25488)
120+++ ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 25489)
121@@ -115,10 +115,10 @@
122 /*Proceed now to heads computations*/
123 solutionsequence_hydro_nonlinear(femmodel);
124 /*If we have a sub-timestep we store the substep inputs in a transient input here*/
125- femmodel->StackTransientInputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput);
126+ femmodel->StackTransientInputonBasex(&substepinput[0],&transientinput[0],subtime,numaveragedinput);
127 }
128 /*averaging the stack*/
129- femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,hydro_averaging);
130+ femmodel->AverageTransientInputonBasex(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,hydro_averaging);
131 /*And reseting to global time*/
132 femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
133 femmodel->parameters->SetParam(global_time,TimeEnum);
134Index: ../trunk-jpl/src/c/classes/FemModel.h
135===================================================================
136--- ../trunk-jpl/src/c/classes/FemModel.h (revision 25488)
137+++ ../trunk-jpl/src/c/classes/FemModel.h (revision 25489)
138@@ -178,7 +178,9 @@
139 void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count);
140 void InitTransientInputx(int* transientinput_enum,int numoutputs);
141 void StackTransientInputx(int* input_enum,int* transientinput_enum,IssmDouble hydrotime,int numoutputs);
142+ void StackTransientInputonBasex(int* input_enum,int* transientinput_enum,IssmDouble hydrotime,int numoutputs);
143 void AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs,int averaging_method);
144+ void AverageTransientInputonBasex(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs,int averaging_method);
145 void UpdateConstraintsx(void);
146 int UpdateVertexPositionsx(void);
147
148Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
149===================================================================
150--- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 25488)
151+++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 25489)
152@@ -687,8 +687,8 @@
153
154 bool element_active;
155 Element* element=NULL;
156-
157- for(int i=0;i<femmodel->elements->Size();i++){
158+ int elementssize=femmodel->elements->Size();
159+ for(int i=0;i<elementssize;i++){
160 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
161
162 Input* input=element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(input);
163@@ -747,8 +747,8 @@
164
165 bool element_active;
166 Element* element=NULL;
167-
168- for(int i=0;i<femmodel->elements->Size();i++){
169+ int elementssize = femmodel->elements->Size();
170+ for(int i=0;i<elementssize;i++){
171 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
172
173 Input* input=element->GetInput(HydrologydcMaskThawedNodeEnum); _assert_(input);
174Index: ../trunk-jpl/src/c/classes/FemModel.cpp
175===================================================================
176--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 25488)
177+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 25489)
178@@ -48,7 +48,6 @@
179 #include "../modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h"
180 #include "../modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h"
181 #include "../modules/NodalValuex/NodalValuex.h"
182-#include "../modules/GetVectorFromInputsx/GetVectorFromInputsx.h"
183 #include "../modules/AverageOntoPartitionx/AverageOntoPartitionx.h"
184 /*}}}*/
185
186@@ -4542,9 +4541,9 @@
187 int npart;
188
189 /*retrieve partition vectors for responses that are scaled:*/
190- this->parameters->FindParam(&response_partitions,&response_partitions_num,NULL,NULL,QmuResponsePartitionsEnum);
191- this->parameters->FindParam(&response_partitions_npart,NULL,NULL,QmuResponsePartitionsNpartEnum);
192-
193+ this->parameters->FindParam(&response_partitions,&response_partitions_num,NULL,NULL,QmuResponsePartitionsEnum);
194+ this->parameters->FindParam(&response_partitions_npart,NULL,NULL,QmuResponsePartitionsNpartEnum);
195+
196 /*retrieve my_rank: */
197 my_rank=IssmComm::GetRank();
198
199@@ -5020,7 +5019,7 @@
200 recurence=new Vector<IssmDouble>(numnodes);
201 this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
202 this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum);
203- GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
204+ GetVectoronBaseFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
205
206 for (int i=0;i<elements->Size();i++){
207 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
208@@ -5067,13 +5066,25 @@
209
210 /*Update node activation accordingly*/
211 int counter = 0; //this is probably not acurate but we are only interested in positivity
212+ Element* basalelement=NULL;
213 for(int i=0;i<elements->Size();i++){
214 Element *element = xDynamicCast<Element*>(elements->GetObjectByOffset(i));
215- int numnodes = element->GetNumberOfNodes();
216+ switch(domaintype){
217+ case Domain2DhorizontalEnum:
218+ basalelement = element;
219+ break;
220+ case Domain3DEnum:
221+ if(!element->IsOnBase()) continue;
222+ basalelement = element->SpawnBasalElement();
223+ break;
224+ default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
225+ }
226+
227+ int numnodes = basalelement->GetNumberOfNodes();
228 IssmDouble *base = xNew<IssmDouble>(numnodes);
229- element->GetInputListOnNodes(&base[0],BaseEnum);
230+ basalelement->GetInputListOnNodes(&base[0],BaseEnum);
231 for(int in=0;in<numnodes;in++){
232- Node* node=element->GetNode(in);
233+ Node* node=basalelement->GetNode(in);
234 if(serial_active[node->Sid()]==1.){
235 node->Activate();
236 if(!node->IsClone()) counter++;
237@@ -5084,6 +5095,10 @@
238 }
239 }
240 xDelete<IssmDouble>(base);
241+ if(domaintype!=Domain2DhorizontalEnum){
242+ basalelement->DeleteMaterials();
243+ delete basalelement;
244+ }
245 }
246 xDelete<IssmDouble>(serial_active);
247 delete effanalysis;
248@@ -5133,7 +5148,7 @@
249 }
250 /*for other cases we just grab the mask from the initialisation value*/
251 else{
252- GetVectorFromInputsx(&serial_mask,this,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
253+ GetVectoronBaseFromInputsx(&serial_mask,this,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
254 }
255 /*Update Mask and elementize*/
256 InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
257@@ -5280,6 +5295,61 @@
258 }
259 }
260 /*}}}*/
261+void FemModel::StackTransientInputonBasex(int* input_enum,int* transientinput_enum,IssmDouble subtime,int numoutputs){ /*{{{*/
262+
263+ Element* basalelement=NULL;
264+ int domaintype;
265+ this->parameters->FindParam(&domaintype,DomainTypeEnum);
266+
267+ for(int i=0;i<numoutputs;i++){
268+ if(input_enum[i]<0){
269+ _error_("Can't deal with non enum fields for result Stack");
270+ }
271+ else{
272+ for(int j=0;j<elements->Size();j++){
273+ /*Get the right transient input*/
274+ Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
275+ /*Get basal element*/
276+ switch(domaintype){
277+ case Domain2DhorizontalEnum:
278+ break;
279+ case Domain3DEnum:
280+ if(!element->IsOnBase()) continue;
281+ break;
282+ default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
283+ }
284+ /*Get values and lid list*/
285+ TransientInput* transientinput = this->inputs->GetTransientInput(transientinput_enum[i]);
286+ const int numvertices = element->GetNumberOfVertices();
287+ IssmDouble* values = xNew<IssmDouble>(numvertices);
288+ int *vertexlids = xNew<int>(numvertices);
289+ switch(domaintype){
290+ case Domain2DhorizontalEnum:
291+ element->GetInputListOnVertices(&values[0],input_enum[i]); //this is the enum to stack
292+ element->GetVerticesLidList(vertexlids);
293+ transientinput->AddTriaTimeInput(subtime,numvertices,vertexlids,values,P1Enum);
294+ break;
295+ case Domain3DEnum:{
296+ element->GetInputListOnVertices(&values[0],input_enum[i]); //this is the enum to stack
297+ element->GetVerticesLidList(vertexlids);
298+ Penta* penta=xDynamicCast<Penta*>(elements->GetObjectByOffset(j));
299+ for(;;){
300+ transientinput->AddPentaTimeInput(subtime,numvertices,vertexlids,values,P1Enum);
301+ if (penta->IsOnSurface()) break;
302+ penta=penta->GetUpperPenta(); _assert_(penta->Id()!=element->id);
303+ penta->GetVerticesLidList(vertexlids);
304+ }
305+ break;
306+ }
307+ default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
308+ }
309+ xDelete<IssmDouble>(values);
310+ xDelete<int>(vertexlids);
311+ }
312+ }
313+ }
314+}
315+/*}}}*/
316 void FemModel::AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs, int averaging_method){ /*{{{*/
317
318 for(int i=0;i<numoutputs;i++){
319@@ -5289,6 +5359,28 @@
320 }
321 }
322 }/*}}}*/
323+void FemModel::AverageTransientInputonBasex(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs, int averaging_method){ /*{{{*/
324+
325+ int domaintype;
326+ this->parameters->FindParam(&domaintype,DomainTypeEnum);
327+
328+ for(int i=0;i<numoutputs;i++){
329+ for(int j=0;j<this->elements->Size();j++){
330+ Element* element = xDynamicCast<Element*>(elements->GetObjectByOffset(j));
331+ /*Get basal element*/
332+ switch(domaintype){
333+ case Domain2DhorizontalEnum:
334+ element->CreateInputTimeAverage(transientinput_enum[i],averagedinput_enum[i],init_time,end_time,averaging_method);;
335+ break;
336+ case Domain3DEnum:
337+ if(!element->IsOnBase()) continue;
338+ element->CreateInputTimeAverage(transientinput_enum[i],averagedinput_enum[i],init_time,end_time,averaging_method);
339+ break;
340+ default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
341+ }
342+ }
343+ }
344+}/*}}}*/
345 #ifdef _HAVE_JAVASCRIPT_
346 FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/
347 /*configuration: */
Note: See TracBrowser for help on using the repository browser.