source: issm/oecreview/Archive/23390-24306/ISSM-23543-23544.diff

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

NEW: adding Archive/23390-24306
File size: 12.2 KB
RevLine 
[24307]1Index: ../trunk-jpl/src/c/cores/hydrology_core.cpp
2===================================================================
3--- ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 23543)
4+++ ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 23544)
5@@ -74,7 +74,8 @@
6 int inputtostack[4]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum,EplHeadHydrostepEnum,HydrologydcEplThicknessHydrostepEnum};
7 int stackedinput[4]={EffectivePressureStackedEnum,SedimentHeadStackedEnum,EplHeadStackedEnum,HydrologydcEplThicknessStackedEnum};
8 int averagedinput[4]={EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum};
9- femmodel->InitTransientOutputx(&stackedinput[0],4);
10+ //femmodel->InitTransientOutputx(&stackedinput[0],4);
11+ femmodel->InitMeanOutputx(&stackedinput[0],4);
12 while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
13 hydrostep+=1;
14 hydrotime+=hydrodt;
15@@ -87,15 +88,18 @@
16 /*Proceed now to heads computations*/
17 solutionsequence_hydro_nonlinear(femmodel);
18 /*If we have a sub-timestep we stack the variables here*/
19- femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,4);
20+ //femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,4);
21+ femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],4);
22 }
23- femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,4);
24+ //femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,4);
25+ femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],4);
26 }
27 else{
28 int inputtostack[2]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum};
29 int stackedinput[2]={EffectivePressureStackedEnum,SedimentHeadStackedEnum};
30 int averagedinput[2]={EffectivePressureEnum,SedimentHeadEnum};
31- femmodel->InitTransientOutputx(&stackedinput[0],2);
32+ //femmodel->InitTransientOutputx(&stackedinput[0],2);
33+ femmodel->InitMeanOutputx(&stackedinput[0],2);
34 while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
35 hydrotime+=hydrodt;
36 /*save preceding timestep*/
37@@ -103,9 +107,11 @@
38 /*Proceed now to heads computations*/
39 solutionsequence_hydro_nonlinear(femmodel);
40 /*If we have a sub-timestep we stack the variables here*/
41- femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,2);
42+ //femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,2);
43+ femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],2);
44 }
45- femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,2);
46+ //femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,2);
47+ femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],2);
48 }
49 }
50 else{
51Index: ../trunk-jpl/src/c/classes/FemModel.cpp
52===================================================================
53--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23543)
54+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23544)
55@@ -4820,15 +4820,15 @@
56 if(VerboseSolution()) _printf0_(" Number of active nodes L2 Projection: "<< counter <<"\n");
57 }
58 /*}}}*/
59-void FemModel::InitTransientOutputx(int* input_enum,int numoutputs){ /*{{{*/
60+void FemModel::InitTransientOutputx(int* stackedinput_enum,int numoutputs){ /*{{{*/
61
62 for(int i=0;i<numoutputs;i++){
63- if(input_enum[i]<0){
64+ if(stackedinput_enum[i]<0){
65 _error_("Can't deal with non enum fields for result Stack");
66 }
67 else{
68 for(int j=0;j<elements->Size();j++){
69- TransientInput* transient_input = new TransientInput(input_enum[i]);
70+ TransientInput* transient_input = new TransientInput(stackedinput_enum[i]);
71 /*Intermediaries*/
72 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
73 //transient_input->Configure(element->parameters);
74@@ -4849,13 +4849,12 @@
75 /*Intermediaries*/
76 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
77 TransientInput* stacking_input=NULL;
78- Input* input=element->inputs->GetInput(stackedinput_enum[i]); _assert_(input);
79- Input* input_to_stack=element->GetInput(input_enum[i]); _assert_(input_to_stack);
80+ Input* input=element->inputs->GetInput(stackedinput_enum[i]); _assert_(input); //this is the enum stack
81 stacking_input=dynamic_cast<TransientInput*>(input);
82
83 int numvertices = element->GetNumberOfVertices();
84 IssmDouble* N=xNew<IssmDouble>(numvertices);
85- element->GetInputListOnVertices(&N[0],input_enum[i]);
86+ element->GetInputListOnVertices(&N[0],input_enum[i]); //this is the enum to stack
87 switch(element->ObjectEnum()){
88 case TriaEnum:
89 stacking_input->AddTimeInput(new TriaInput(stackedinput_enum[i],&N[0],P1Enum),hydrotime);
90@@ -4874,10 +4873,10 @@
91 }
92 }
93 /*}}}*/
94-void FemModel::AverageTransientOutputx(int* input_enum,int* averagedinput_enum,IssmDouble init_time,int numoutputs){ /*{{{*/
95+void FemModel::AverageTransientOutputx(int* stackedinput_enum,int* averagedinput_enum,IssmDouble init_time,int numoutputs){ /*{{{*/
96
97 for(int i=0;i<numoutputs;i++){
98- if(input_enum[i]<0){
99+ if(stackedinput_enum[i]<0){
100 _error_("Can't deal with non enum fields for result Stack");
101 }
102 else{
103@@ -4887,7 +4886,7 @@
104 int numvertices = element->GetNumberOfVertices();
105 IssmDouble* time_averaged=NULL;
106
107- Input* input=element->inputs->GetInput(input_enum[i]); _assert_(input);
108+ Input* input=element->inputs->GetInput(stackedinput_enum[i]); _assert_(input);
109 TransientInput* stacking_input=NULL;
110 stacking_input=dynamic_cast<TransientInput*>(input);
111 stacking_input->GetInputAverageOnTimes(&time_averaged,init_time);
112@@ -4899,6 +4898,95 @@
113 }
114 }
115 /*}}}*/
116+void FemModel::InitMeanOutputx(int* stackedinput_enum,int numoutputs){ /*{{{*/
117+
118+ for(int i=0;i<numoutputs;i++){
119+ if(stackedinput_enum[i]<0){
120+ _error_("Can't deal with non enum fields for result Stack");
121+ }
122+ else{
123+ for(int j=0;j<elements->Size();j++){
124+ /*Intermediaries*/
125+ Element* element =xDynamicCast<Element*>(elements->GetObjectByOffset(j));
126+ int numvertices =element->GetNumberOfVertices();
127+ IssmDouble zeros[numvertices] ={0.0};
128+ switch(element->ObjectEnum()){
129+ case TriaEnum:
130+ element->inputs->AddInput(new TriaInput(stackedinput_enum[i],&zeros[0],P1Enum));
131+ break;
132+ case PentaEnum:
133+ element->inputs->AddInput(new PentaInput(stackedinput_enum[i],&zeros[0],P1Enum));
134+ break;
135+ case TetraEnum:
136+ element->inputs->AddInput(new TetraInput(stackedinput_enum[i],&zeros[0],P1Enum));
137+ break;
138+ default: _error_("Not implemented yet");
139+ }
140+ }
141+ }
142+ }
143+}
144+/*}}}*/
145+void FemModel::SumOutputx(int* input_enum,int* stackedinput_enum,int numoutputs){ /*{{{*/
146+
147+ //First get sub-timestep
148+ IssmDouble hydrodt;
149+ this->parameters->FindParam(&hydrodt,HydrologydtEnum);
150+
151+ for(int i=0;i<numoutputs;i++){
152+ if(input_enum[i]<0){
153+ _error_("Can't deal with non enum fields for result Stack");
154+ }
155+ else{
156+ for(int j=0;j<elements->Size();j++){
157+ /*Intermediaries*/
158+ Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
159+ int numvertices = element->GetNumberOfVertices();
160+ IssmDouble* values_to_add=xNew<IssmDouble>(numvertices);
161+ IssmDouble* existing_values=xNew<IssmDouble>(numvertices);
162+ element->GetInputListOnVertices(&values_to_add[0],input_enum[i]); //those are the values to add
163+ element->GetInputListOnVertices(&existing_values[0],stackedinput_enum[i]); //those are the values to add
164+ for(int k=0;k<numvertices;k++){
165+ existing_values[k]+=values_to_add[k]*hydrodt;
166+ }
167+ element->AddInput(stackedinput_enum[i],&existing_values[0],P1Enum);
168+ xDelete<IssmDouble>(existing_values);
169+ xDelete<IssmDouble>(values_to_add);
170+ }
171+ }
172+ }
173+}
174+/*}}}*/
175+void FemModel::AverageSumOutputx(int* stackedinput_enum,int* averagedinput_enum,int numoutputs){ /*{{{*/
176+
177+ //First get timestep
178+ IssmDouble maindt;
179+ this->parameters->FindParam(&maindt,TimesteppingTimeStepEnum);
180+ for(int i=0;i<numoutputs;i++){
181+ if(stackedinput_enum[i]<0){
182+ _error_("Can't deal with non enum fields for result Stack");
183+ }
184+ else{
185+ for(int j=0;j<elements->Size();j++){
186+ /*Intermediaries*/
187+ Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
188+ int numvertices = element->GetNumberOfVertices();
189+ IssmDouble* time_averaged=xNew<IssmDouble>(numvertices);
190+ IssmDouble* existing_values=xNew<IssmDouble>(numvertices);
191+ element->GetInputListOnVertices(&existing_values[0],stackedinput_enum[i]); //those are the values to add
192+
193+ for(int k=0;k<numvertices;k++){
194+ time_averaged[k]=existing_values[k]/maindt;
195+ }
196+
197+ element->AddInput(averagedinput_enum[i],&time_averaged[0],P1Enum);
198+ xDelete<IssmDouble>(time_averaged);
199+ xDelete<IssmDouble>(existing_values);
200+ }
201+ }
202+ }
203+}
204+/*}}}*/
205 #ifdef _HAVE_JAVASCRIPT_
206 FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/
207 /*configuration: */
208Index: ../trunk-jpl/src/c/classes/FemModel.h
209===================================================================
210--- ../trunk-jpl/src/c/classes/FemModel.h (revision 23543)
211+++ ../trunk-jpl/src/c/classes/FemModel.h (revision 23544)
212@@ -140,14 +140,14 @@
213 void Deflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y);
214 #endif
215 #ifdef _HAVE_ESA_
216- void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy);
217- void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
218+ void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy);
219+ void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
220 #endif
221 #ifdef _HAVE_SEALEVELRISE_
222 void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop);
223 void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
224 void SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
225- void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
226+ void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
227 IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg);
228 #endif
229 void HydrologyEPLupdateDomainx(IssmDouble* pEplcount);
230@@ -158,7 +158,10 @@
231 void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count);
232 void InitTransientOutputx(int* input_enum, int numoutputs);
233 void StackTransientOutputx(int* input_enum,int* stackedinput_enum, IssmDouble hydrotime, int numoutputs);
234- void AverageTransientOutputx(int* input_enum,int* averagedinput_enum,IssmDouble hydrotime,int numoutputs);
235+ void AverageTransientOutputx(int* stackedinput_enum,int* averagedinput_enum,IssmDouble init_time,int numoutputs);
236+ void InitMeanOutputx(int* stackedinput_enum,int numoutputs);
237+ void SumOutputx(int* input_enum,int* stackedinput_enum,int numoutputs);
238+ void AverageSumOutputx(int* stackedinput_enum,int* averagedinput_enum,int numoutputs);
239 void UpdateConstraintsx(void);
240 int UpdateVertexPositionsx(void);
241
Note: See TracBrowser for help on using the repository browser.