source: issm/oecreview/Archive/23390-24306/ISSM-23887-23888.diff@ 24307

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

NEW: adding Archive/23390-24306
File size: 15.2 KB
RevLine 
[24307]1Index: ../trunk-jpl/src/c/classes/FemModel.cpp
2===================================================================
3--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23887)
4+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23888)
5@@ -1545,6 +1545,22 @@
6 *pM=total_mass_flux;
7
8 }/*}}}*/
9+void FemModel::GroundinglineMassFluxx(IssmDouble* pM, bool scaled){/*{{{*/
10+
11+ IssmDouble local_mass_flux = 0;
12+ IssmDouble total_mass_flux;
13+
14+ for(int i=0;i<this->elements->Size();i++){
15+ Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
16+ local_mass_flux+=element->GroundinglineMassFlux(scaled);
17+ }
18+ ISSM_MPI_Reduce(&local_mass_flux,&total_mass_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
19+ ISSM_MPI_Bcast(&total_mass_flux,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
20+
21+ /*Assign output pointers: */
22+ *pM=total_mass_flux;
23+
24+}/*}}}*/
25 void FemModel::IceMassx(IssmDouble* pM, bool scaled){/*{{{*/
26
27 IssmDouble local_ice_mass = 0;
28@@ -2205,6 +2221,7 @@
29 case IceVolumeAboveFloatationScaledEnum: this->IceVolumeAboveFloatationx(&double_result,true); break;
30 case GroundedAreaEnum: this->GroundedAreax(&double_result,false); break;
31 case GroundedAreaScaledEnum: this->GroundedAreax(&double_result,true); break;
32+ case GroundinglineMassFluxEnum: this->GroundinglineMassFluxx(&double_result,false); break;
33 case FloatingAreaEnum: this->FloatingAreax(&double_result,false); break;
34 case FloatingAreaScaledEnum: this->FloatingAreax(&double_result,true); break;
35 case MinVelEnum: this->MinVelx(&double_result); break;
36@@ -2402,7 +2419,7 @@
37 case IceVolumeScaledEnum: this->IceVolumex(responses, true); break;
38 case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(responses, false); break;
39 case IceVolumeAboveFloatationScaledEnum: this->IceVolumeAboveFloatationx(responses, true); break;
40- case IcefrontMassFluxEnum: this->IcefrontMassFluxx(responses, true); break;
41+ case IcefrontMassFluxEnum: this->IcefrontMassFluxx(responses, false); break;
42 case GroundedAreaEnum: this->GroundedAreax(responses, false); break;
43 case GroundedAreaScaledEnum: this->GroundedAreax(responses, true); break;
44 case FloatingAreaEnum: this->FloatingAreax(responses, false); break;
45Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
46===================================================================
47--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 23887)
48+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 23888)
49@@ -94,6 +94,8 @@
50 bool HasEdgeOnSurface();
51 IssmDouble IceVolume(bool scaled);
52 IssmDouble IceVolumeAboveFloatation(bool scaled);
53+ IssmDouble IcefrontMassFlux(bool scaled);
54+ IssmDouble GroundinglineMassFlux(bool scaled);
55 void Ismip6FloatingiceMeltingRate(void);
56 void InputDepthAverageAtBase(int enum_type,int average_enum_type);
57 void InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/};
58Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
59===================================================================
60--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23887)
61+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23888)
62@@ -1079,6 +1079,9 @@
63 IssmDouble* H=NULL;;
64 int nrfrontbed,numiceverts;
65
66+ if(!this->IsOnBase()) return 0;
67+ if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
68+
69 /*Retrieve all inputs and parameters*/
70 GetInputListOnVertices(&bed[0],BedEnum);
71 GetInputListOnVertices(&surfaces[0],SurfaceEnum);
72@@ -1085,9 +1088,6 @@
73 GetInputListOnVertices(&bases[0],BaseEnum);
74 GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
75
76- if(!this->IsOnBase()) return 0;
77- if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
78-
79 nrfrontbed=0;
80 for(int i=0;i<NUMVERTICES2D;i++){
81 /*Find if bed<0*/
82Index: ../trunk-jpl/src/c/classes/Elements/Element.h
83===================================================================
84--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 23887)
85+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 23888)
86@@ -240,6 +240,7 @@
87 virtual IssmDouble IceVolume(bool scaled)=0;
88 virtual IssmDouble IceVolumeAboveFloatation(bool scaled)=0;
89 virtual IssmDouble IcefrontMassFlux(bool scaled){_error_("not implemented");};
90+ virtual IssmDouble GroundinglineMassFlux(bool scaled){_error_("not implemented");};
91 virtual void InputDepthAverageAtBase(int enum_type,int average_enum_type)=0;
92 virtual void InputExtrude(int input_enum,int start)=0;
93 virtual void InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum)=0;
94Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
95===================================================================
96--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23887)
97+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23888)
98@@ -1408,6 +1408,8 @@
99 IssmDouble* H=NULL;;
100 int nrfrontbed,numiceverts;
101
102+ if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
103+
104 /*Retrieve all inputs and parameters*/
105 GetInputListOnVertices(&bed[0],BedEnum);
106 GetInputListOnVertices(&surfaces[0],SurfaceEnum);
107@@ -1414,8 +1416,6 @@
108 GetInputListOnVertices(&bases[0],BaseEnum);
109 GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
110
111- if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
112-
113 nrfrontbed=0;
114 for(int i=0;i<NUMVERTICES;i++){
115 /*Find if bed<0*/
116@@ -1946,6 +1946,262 @@
117 }
118 }
119 /*}}}*/
120+IssmDouble Tria::IcefrontMassFlux(bool scaled){/*{{{*/
121+
122+ /*Make sure there is an ice front here*/
123+ if(!IsIceInElement() || !IsIcefront()) return 0;
124+
125+ /*Scaled not implemented yet...*/
126+ _assert_(!scaled);
127+
128+ int domaintype,index1,index2;
129+ const IssmPDouble epsilon = 1.e-15;
130+ IssmDouble s1,s2;
131+ IssmDouble gl[NUMVERTICES];
132+ IssmDouble xyz_front[2][3];
133+
134+
135+ IssmDouble *xyz_list = NULL;
136+ this->GetVerticesCoordinates(&xyz_list);
137+
138+ /*Recover parameters and values*/
139+ parameters->FindParam(&domaintype,DomainTypeEnum);
140+ GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
141+
142+ /*Be sure that values are not zero*/
143+ if(gl[0]==0.) gl[0]=gl[0]+epsilon;
144+ if(gl[1]==0.) gl[1]=gl[1]+epsilon;
145+ if(gl[2]==0.) gl[2]=gl[2]+epsilon;
146+
147+ if(domaintype==Domain2DverticalEnum){
148+ _error_("not implemented");
149+ }
150+ else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
151+ int pt1 = 0;
152+ int pt2 = 1;
153+ if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
154+
155+ /*Portion of the segments*/
156+ s1=gl[2]/(gl[2]-gl[1]);
157+ s2=gl[2]/(gl[2]-gl[0]);
158+ if(gl[2]<0.){
159+ pt1 = 1; pt2 = 0;
160+ }
161+ xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
162+ xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
163+ xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
164+ xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
165+ xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
166+ xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
167+ }
168+ else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
169+
170+ /*Portion of the segments*/
171+ s1=gl[0]/(gl[0]-gl[1]);
172+ s2=gl[0]/(gl[0]-gl[2]);
173+ if(gl[0]<0.){
174+ pt1 = 1; pt2 = 0;
175+ }
176+
177+ xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
178+ xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
179+ xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
180+ xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
181+ xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
182+ xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
183+ }
184+ else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
185+
186+ /*Portion of the segments*/
187+ s1=gl[1]/(gl[1]-gl[0]);
188+ s2=gl[1]/(gl[1]-gl[2]);
189+ if(gl[1]<0.){
190+ pt1 = 1; pt2 = 0;
191+ }
192+
193+ xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
194+ xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
195+ xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
196+ xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
197+ xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
198+ xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
199+ }
200+ else{
201+ _error_("case not possible");
202+ }
203+
204+ }
205+ else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
206+
207+ /*Some checks in debugging mode*/
208+ _assert_(s1>=0 && s1<=1.);
209+ _assert_(s2>=0 && s2<=1.);
210+
211+ /*Get normal vector*/
212+ IssmDouble normal[3];
213+ this->NormalSection(&normal[0],&xyz_front[0][0]);
214+ normal[0] = -normal[0];
215+ normal[1] = -normal[1];
216+
217+ /*Get inputs*/
218+ IssmDouble flux = 0.;
219+ IssmDouble vx,vy,thickness,Jdet;
220+ IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
221+ Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
222+ Input* vx_input=NULL;
223+ Input* vy_input=NULL;
224+ if(domaintype==Domain2DhorizontalEnum){
225+ vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
226+ vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
227+ }
228+ else{
229+ vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
230+ vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
231+ }
232+
233+ /*Start looping on Gaussian points*/
234+ Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
235+ for(int ig=gauss->begin();ig<gauss->end();ig++){
236+
237+ gauss->GaussPoint(ig);
238+ thickness_input->GetInputValue(&thickness,gauss);
239+ vx_input->GetInputValue(&vx,gauss);
240+ vy_input->GetInputValue(&vy,gauss);
241+ this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
242+
243+ flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
244+ }
245+
246+ return flux;
247+}
248+/*}}}*/
249+IssmDouble Tria::GroundinglineMassFlux(bool scaled){/*{{{*/
250+
251+ /*Make sure there is an ice front here*/
252+ if(!IsIceInElement()) return 0;
253+ if(!IsZeroLevelset(MaskGroundediceLevelsetEnum)) return 0;
254+
255+ /*Scaled not implemented yet...*/
256+ _assert_(!scaled);
257+
258+ int domaintype,index1,index2;
259+ const IssmPDouble epsilon = 1.e-15;
260+ IssmDouble s1,s2;
261+ IssmDouble gl[NUMVERTICES];
262+ IssmDouble xyz_front[2][3];
263+
264+ IssmDouble *xyz_list = NULL;
265+ this->GetVerticesCoordinates(&xyz_list);
266+
267+ /*Recover parameters and values*/
268+ parameters->FindParam(&domaintype,DomainTypeEnum);
269+ GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
270+
271+ /*Be sure that values are not zero*/
272+ if(gl[0]==0.) gl[0]=gl[0]+epsilon;
273+ if(gl[1]==0.) gl[1]=gl[1]+epsilon;
274+ if(gl[2]==0.) gl[2]=gl[2]+epsilon;
275+
276+ if(domaintype==Domain2DverticalEnum){
277+ _error_("not implemented");
278+ }
279+ else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
280+ int pt1 = 0;
281+ int pt2 = 1;
282+ if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
283+
284+ /*Portion of the segments*/
285+ s1=gl[2]/(gl[2]-gl[1]);
286+ s2=gl[2]/(gl[2]-gl[0]);
287+ if(gl[2]<0.){
288+ pt1 = 1; pt2 = 0;
289+ }
290+ xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
291+ xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
292+ xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
293+ xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
294+ xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
295+ xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
296+ }
297+ else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
298+
299+ /*Portion of the segments*/
300+ s1=gl[0]/(gl[0]-gl[1]);
301+ s2=gl[0]/(gl[0]-gl[2]);
302+ if(gl[0]<0.){
303+ pt1 = 1; pt2 = 0;
304+ }
305+
306+ xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
307+ xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
308+ xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
309+ xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
310+ xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
311+ xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
312+ }
313+ else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
314+
315+ /*Portion of the segments*/
316+ s1=gl[1]/(gl[1]-gl[0]);
317+ s2=gl[1]/(gl[1]-gl[2]);
318+ if(gl[1]<0.){
319+ pt1 = 1; pt2 = 0;
320+ }
321+
322+ xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
323+ xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
324+ xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
325+ xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
326+ xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
327+ xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
328+ }
329+ else{
330+ _error_("case not possible");
331+ }
332+
333+ }
334+ else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
335+
336+ /*Some checks in debugging mode*/
337+ _assert_(s1>=0 && s1<=1.);
338+ _assert_(s2>=0 && s2<=1.);
339+
340+ /*Get normal vector*/
341+ IssmDouble normal[3];
342+ this->NormalSection(&normal[0],&xyz_front[0][0]);
343+
344+ /*Get inputs*/
345+ IssmDouble flux = 0.;
346+ IssmDouble vx,vy,thickness,Jdet;
347+ IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
348+ Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
349+ Input* vx_input=NULL;
350+ Input* vy_input=NULL;
351+ if(domaintype==Domain2DhorizontalEnum){
352+ vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
353+ vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
354+ }
355+ else{
356+ vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
357+ vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
358+ }
359+
360+ /*Start looping on Gaussian points*/
361+ Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
362+ for(int ig=gauss->begin();ig<gauss->end();ig++){
363+
364+ gauss->GaussPoint(ig);
365+ thickness_input->GetInputValue(&thickness,gauss);
366+ vx_input->GetInputValue(&vx,gauss);
367+ vy_input->GetInputValue(&vy,gauss);
368+ this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
369+
370+ flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
371+ }
372+
373+ return flux;
374+}
375+/*}}}*/
376 IssmDouble Tria::IceVolume(bool scaled){/*{{{*/
377
378 /*The volume of a truncated prism is area_base * 1/numedges sum(length of edges)*/
379Index: ../trunk-jpl/src/c/classes/FemModel.h
380===================================================================
381--- ../trunk-jpl/src/c/classes/FemModel.h (revision 23887)
382+++ ../trunk-jpl/src/c/classes/FemModel.h (revision 23888)
383@@ -101,6 +101,7 @@
384 void GroundedAreax(IssmDouble* pV, bool scaled);
385 void IcefrontAreax();
386 void IcefrontMassFluxx(IssmDouble* presponse, bool scaled);
387+ void GroundinglineMassFluxx(IssmDouble* presponse, bool scaled);
388 void IceMassx(IssmDouble* pV, bool scaled);
389 void IceVolumex(IssmDouble* pV, bool scaled);
390 void IceVolumeAboveFloatationx(IssmDouble* pV, bool scaled);
Note: See TracBrowser for help on using the repository browser.