source: issm/oecreview/Archive/25834-26739/ISSM-26652-26653.diff@ 26740

Last change on this file since 26740 was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739
File size: 9.4 KB
RevLine 
[26740]1Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
2===================================================================
3--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 26652)
4+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 26653)
5@@ -1295,7 +1295,7 @@
6 else phi = f1*f2;
7
8 /*Compute weights*/
9- gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2);
10+ gauss = this->NewGauss(point,f1,f2,1-trapezeisnegative,2);
11
12 total_weight = 0.0;
13 for(int i=0;i<NUMVERTICES2D;i++)weights[i] = 0;
14@@ -1305,10 +1305,8 @@
15 total_weight += gauss->weight;
16 }
17
18- /*Normalize to phi*/
19- if(total_weight>0.){
20- for(int i=0;i<NUMVERTICES2D;i++) weights[i] = weights[i]*phi/total_weight;
21- }
22+ /*Normalizing to phi such that weights provide coefficients for integration over subelement (for averaging:phi*weights)*/
23+ if(total_weight>0.) for(int i=0;i<NUMVERTICES;i++) weights[i] = weights[i]*phi/total_weight;
24 else for(int i=0;i<NUMVERTICES2D;i++) weights[i] = 0.0;
25
26 /*Cleanup*/
27@@ -2206,7 +2204,7 @@
28
29 IssmDouble basetot;
30 height = 0.0;
31- for(int i=0;i<NUMVERTICES2D;i++) height += weights[i]*heights[i];
32+ for(int i=0;i<NUMVERTICES2D;i++) height += weights[i]/phi*heights[i];
33 basetot = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
34 base = basetot*phi;
35
36@@ -2216,7 +2214,7 @@
37 Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
38 /*Compute loop only over lower vertices: i<NUMVERTICES2D*/
39 scalefactor = 0.0;
40- for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]*scalefactor_vertices[i];
41+ for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
42 base = base*scalefactor;
43 xDelete<IssmDouble>(scalefactor_vertices);
44 }
45@@ -4388,6 +4386,7 @@
46 /*The smb[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
47 IssmDouble base,smb,rho_ice,scalefactor;
48 IssmDouble Total_Smb=0;
49+ IssmDouble lsf[NUMVERTICES];
50 IssmDouble xyz_list[NUMVERTICES][3];
51
52 /*Get material parameters :*/
53@@ -4403,16 +4402,47 @@
54 base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
55
56 /*Now get the average SMB over the element*/
57- Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
58+ Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
59+ if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
60+ /*Partially ice-covered element*/
61+ bool mainlyice;
62+ int point;
63+ IssmDouble* smb_vertices = xNew<IssmDouble>(NUMVERTICES);
64+ IssmDouble weights[NUMVERTICES2D];
65+ IssmDouble lsf2d[NUMVERTICES2D];
66+ IssmDouble f1,f2,phi;
67+ Element::GetInputListOnVertices(&smb_vertices[0],SmbMassBalanceEnum);
68+ for(int i=0;i<NUMVERTICES2D;i++) lsf2d[i] = lsf[i];
69+ GetFractionGeometry2D(weights,&phi,&point,&f1,&f2,&mainlyice,lsf2d);
70+ smb = 0.0;
71+ for(int i=0;i<NUMVERTICES2D;i++) smb += weights[i]*smb_vertices[i];
72
73- smb_input->GetInputAverage(&smb);
74- if(scaled==true){
75- Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
76- scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
77+ if(scaled==true){
78+ IssmDouble* scalefactor_vertices = xNew<IssmDouble>(NUMVERTICES);
79+ Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
80+ /*Compute loop only over lower vertices: i<NUMVERTICES2D*/
81+ scalefactor = 0.0;
82+ for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
83+ xDelete<IssmDouble>(scalefactor_vertices);
84+ }
85+ else scalefactor = 1.0;
86+
87+ /*Cleanup*/
88+ xDelete<IssmDouble>(smb_vertices);
89 }
90+
91 else{
92- scalefactor=1.;
93+ /*Fully ice-covered element*/
94+ Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
95+ smb_input->GetInputAverage(&smb);
96+
97+ if(scaled==true){
98+ Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
99+ scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
100+ }
101+ else scalefactor=1.0;
102 }
103+
104 Total_Smb=rho_ice*base*smb*scalefactor;// smb on element in kg s-1
105
106 /*Return: */
107Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
108===================================================================
109--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26652)
110+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26653)
111@@ -1816,7 +1816,7 @@
112 else phi=f1*f2;
113
114 /*Compute weights:*/
115- gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2);
116+ gauss = this->NewGauss(point,f1,f2,1-trapezeisnegative,2); //VV correction (16Nov2021)
117
118 total_weight=0;
119 for(int i=0;i<NUMVERTICES;i++)weights[i]=0;
120@@ -1825,7 +1825,7 @@
121 for(int i=0;i<NUMVERTICES;i++)weights[i]+=loadweights_g[i]*gauss->weight;
122 total_weight+=gauss->weight;
123 }
124- //normalize to phi.
125+ /*Normalizing to phi such that weights provide coefficients for integration over subelement (for averaging:phi*weights)*/
126 if(total_weight>0.) for(int i=0;i<NUMVERTICES;i++)weights[i]/=total_weight/phi;
127 else for(int i=0;i<NUMVERTICES;i++)weights[i]=0;
128
129@@ -3545,8 +3545,8 @@
130 Haverage/=IssmDouble(numthk);
131 }*/
132
133- IssmDouble* lsf = xNew<IssmDouble>(NUMVERTICES);
134- Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
135+ IssmDouble lsf[NUMVERTICES];
136+ Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
137 /*Deal with partially ice-covered elements*/
138 if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
139 bool istrapneg;
140@@ -3562,7 +3562,8 @@
141 GetFractionGeometry(weights,&phi,&point,&f1,&f2,&istrapneg,lsf);
142 for(int i=0;i<NUMVERTICES;i++) Hice[i] = surfaces[i]-bases[i];
143 Haverage = 0.0;
144- for(int i=0;i<NUMVERTICES;i++) Haverage += weights[i]*Hice[i];
145+ /*Use weights[i]/phi to get average thickness over subelement*/
146+ for(int i=0;i<NUMVERTICES;i++) Haverage += weights[i]/phi*Hice[i];
147 /*Get back area of ice-covered base*/
148 area_basetot = this->GetArea();
149 area_base = phi*area_basetot;
150@@ -3572,7 +3573,7 @@
151 IssmDouble* scalefactor_vertices = xNew<IssmDouble>(NUMVERTICES);
152 Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
153 scalefactor = 0.0;
154- for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]*scalefactor_vertices[i];
155+ for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
156 area_base = area_base*scalefactor;
157 xDelete<IssmDouble>(scalefactor_vertices);
158 }
159@@ -3605,7 +3606,6 @@
160 xDelete<int>(indices);
161 xDelete<IssmDouble>(H);
162 xDelete<IssmDouble>(SF);
163- xDelete<IssmDouble>(lsf);
164
165 if(domaintype==Domain2DverticalEnum){
166 return area_base;
167@@ -5478,6 +5478,7 @@
168 /*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/
169 IssmDouble base,smb,rho_ice,scalefactor;
170 IssmDouble Total_Smb=0;
171+ IssmDouble lsf[NUMVERTICES];
172 IssmDouble xyz_list[NUMVERTICES][3];
173
174 /*Get material parameters :*/
175@@ -5491,17 +5492,47 @@
176 * http://en.wikipedia.org/wiki/Triangle
177 * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
178 base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); // area of element in m2
179+
180+ /*Now get the average SMB over the element*/
181+ Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
182+ if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
183+ /*Partially ice-covered element*/
184+ bool mainlyice;
185+ int point;
186+ IssmDouble* weights = xNew<IssmDouble>(NUMVERTICES);
187+ IssmDouble* smb_vertices = xNew<IssmDouble>(NUMVERTICES);
188+ IssmDouble f1,f2,phi;
189+
190+ Element::GetInputListOnVertices(&smb_vertices[0],SmbMassBalanceEnum);
191+ GetFractionGeometry(weights,&phi,&point,&f1,&f2,&mainlyice,lsf);
192+ smb = 0.0;
193+ for(int i=0;i<NUMVERTICES;i++) smb += weights[i]*smb_vertices[i];
194+
195+ if(scaled==true){
196+ IssmDouble* scalefactor_vertices = xNew<IssmDouble>(NUMVERTICES);
197+ Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
198+ scalefactor = 0.0;
199+ for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
200+ xDelete<IssmDouble>(scalefactor_vertices);
201+ }
202+ else scalefactor = 1.0;
203
204- /*Now get the average SMB over the element*/
205- Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
206- smb_input->GetInputAverage(&smb); // average smb on element in m ice s-1
207- if(scaled==true){
208- Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
209- scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
210+ /*Cleanup*/
211+ xDelete<IssmDouble>(weights);
212+ xDelete<IssmDouble>(smb_vertices);
213 }
214 else{
215- scalefactor=1.;
216+ /*Fully ice-covered element*/
217+ Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
218+ smb_input->GetInputAverage(&smb); // average smb on element in m ice s-1
219+
220+ if(scaled==true){
221+ Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
222+ scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
223+ }
224+ else scalefactor=1.0;
225 }
226+
227 Total_Smb=rho_ice*base*smb*scalefactor; // smb on element in kg s-1
228
229 /*Return: */
Note: See TracBrowser for help on using the repository browser.