source: issm/oecreview/Archive/23390-24306/ISSM-23800-23801.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: 12.6 KB
RevLine 
[24307]1Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
2===================================================================
3--- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 23800)
4+++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 23801)
5@@ -74,10 +74,12 @@
6 int GetElementType(void);
7 void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
8 IssmDouble GetGroundedPortion(IssmDouble* xyz_list);
9+ IssmDouble GetIcefrontArea();
10 void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
11 void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
12 void GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented yet");};
13 void GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues){_error_("not implemented yet");};
14+ void GetLevelsetIntersectionBase(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level);
15 Node* GetNode(int node_number);
16 int GetNodeIndex(Node* node);
17 int GetNumberOfNodes(void);
18@@ -149,6 +151,7 @@
19 void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
20 void ResetFSBasalBoundaryCondition(void);
21 void ResetHooks();
22+ void RignotMeltParameterization();
23 void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int N,int M);
24 void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
25 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
26Index: ../trunk-jpl/src/c/classes/FemModel.cpp
27===================================================================
28--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23800)
29+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23801)
30@@ -1497,9 +1497,7 @@
31 }/*}}}*/
32 void FemModel::IcefrontAreax(){/*{{{*/
33
34- int numvertices = this->GetElementsWidth();
35 int numbasins;
36- IssmDouble* BasinId = xNew<IssmDouble>(numvertices);
37 this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
38 IssmDouble* basin_icefront_area = xNewZeroInit<IssmDouble>(numbasins);
39
40@@ -1509,24 +1507,27 @@
41
42 for(int i=0;i<this->elements->Size();i++){
43 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
44+ if(!element->IsOnBase()) continue;
45+ int numvertices = element->GetNumberOfVertices();
46+ IssmDouble* BasinId = xNew<IssmDouble>(numvertices);
47 element->GetInputListOnVertices(BasinId,FrontalForcingsBasinIdEnum);
48- for(int j=0;j<numvertices;j++){
49+ for(int j=0;j<3;j++){
50 if(BasinId[j]==basin){
51 local_icefront_area+=element->GetIcefrontArea();
52 break;
53 }
54 }
55+ xDelete<IssmDouble>(BasinId);
56 }
57 ISSM_MPI_Reduce(&local_icefront_area,&total_icefront_area,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
58 ISSM_MPI_Bcast(&total_icefront_area,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
59-
60+
61 basin_icefront_area[basin-1]=total_icefront_area;
62 }
63
64 this->parameters->AddObject(new DoubleVecParam(FrontalForcingsBasinIcefrontAreaEnum,basin_icefront_area,numbasins));
65+ xDelete<IssmDouble>(basin_icefront_area);
66
67- xDelete<IssmDouble>(basin_icefront_area);
68- xDelete<IssmDouble>(BasinId);
69 }/*}}}*/
70 void FemModel::IceMassx(IssmDouble* pM, bool scaled){/*{{{*/
71
72Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
73===================================================================
74--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23800)
75+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23801)
76@@ -1069,6 +1069,107 @@
77 return phi;
78 }
79 /*}}}*/
80+IssmDouble Penta::GetIcefrontArea(){/*{{{*/
81+
82+ IssmDouble bed[NUMVERTICES]; //basinId[NUMVERTICES];
83+ IssmDouble Haverage,frontarea;
84+ IssmDouble x1,y1,x2,y2,distance;
85+ IssmDouble lsf[NUMVERTICES], Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
86+ int* indices=NULL;
87+ IssmDouble* H=NULL;;
88+ int nrfrontbed,numiceverts;
89+
90+ /*Retrieve all inputs and parameters*/
91+ GetInputListOnVertices(&bed[0],BedEnum);
92+ GetInputListOnVertices(&surfaces[0],SurfaceEnum);
93+ GetInputListOnVertices(&bases[0],BaseEnum);
94+ GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
95+
96+ if(!this->IsOnBase()) return 0;
97+ if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
98+
99+ nrfrontbed=0;
100+ for(int i=0;i<NUMVERTICES2D;i++){
101+ /*Find if bed<0*/
102+ if(bed[i]<0.) nrfrontbed++;
103+ }
104+
105+ if(nrfrontbed==3){
106+ /*2. Find coordinates of where levelset crosses 0*/
107+ int numiceverts;
108+ IssmDouble s[2],x[2],y[2];
109+ this->GetLevelsetIntersectionBase(&indices, &numiceverts,&s[0],MaskIceLevelsetEnum,0.);
110+ _assert_(numiceverts);
111+
112+ /*3 Write coordinates*/
113+ IssmDouble xyz_list[NUMVERTICES][3];
114+ ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
115+ int counter = 0;
116+ if((numiceverts>0) && (numiceverts<NUMVERTICES2D)){
117+ for(int i=0;i<numiceverts;i++){
118+ for(int n=numiceverts;n<NUMVERTICES2D;n++){ // iterate over no-ice vertices
119+ x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
120+ y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
121+ counter++;
122+ }
123+ }
124+ }
125+ else if(numiceverts==NUMVERTICES2D){ //NUMVERTICES ice vertices: calving front lies on element edge
126+
127+ for(int i=0;i<NUMVERTICES2D;i++){
128+ if(lsf[indices[i]]==0.){
129+ x[counter]=xyz_list[indices[i]][0];
130+ y[counter]=xyz_list[indices[i]][1];
131+ counter++;
132+ }
133+ if(counter==2) break;
134+ }
135+ if(counter==1){
136+ /*We actually have only 1 vertex on levelset, write a single point as a segment*/
137+ x[counter]=x[0];
138+ y[counter]=y[0];
139+ counter++;
140+ }
141+ }
142+ else{
143+ _error_("not sure what's going on here...");
144+ }
145+ x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
146+ distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));
147+
148+ int numthk=numiceverts+2;
149+ H=xNew<IssmDouble>(numthk);
150+ for(int iv=0;iv<NUMVERTICES2D;iv++) Haux[iv]=-bed[indices[iv]]; //sort bed in ice/noice
151+
152+ switch(numiceverts){
153+ case 1: // average over triangle
154+ H[0]=Haux[0];
155+ H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
156+ H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
157+ Haverage=(H[1]+H[2])/2;
158+ break;
159+ case 2: // average over quadrangle
160+ H[0]=Haux[0];
161+ H[1]=Haux[1];
162+ H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
163+ H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
164+ Haverage=(H[2]+H[3])/2;
165+ break;
166+ default:
167+ _error_("Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
168+ break;
169+ }
170+ frontarea=distance*Haverage;
171+ }
172+ else return 0;
173+
174+ xDelete<int>(indices);
175+ xDelete<IssmDouble>(H);
176+
177+ _assert_(frontarea>0);
178+ return frontarea;
179+}
180+/*}}}*/
181 void Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
182
183 /* Intermediaries */
184@@ -1132,6 +1233,82 @@
185 return lower_penta;
186 }
187 /*}}}*/
188+void Penta::GetLevelsetIntersectionBase(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level){/*{{{*/
189+
190+ /* GetLevelsetIntersection computes:
191+ * 1. indices of element, sorted in [iceverts, noiceverts] in counterclockwise fashion,
192+ * 2. fraction of intersected triangle edges intersected by levelset, lying below level*/
193+
194+ /*Intermediaries*/
195+ int i, numiceverts, numnoiceverts;
196+ int ind0, ind1, lastindex;
197+ int indices_ice[NUMVERTICES2D],indices_noice[NUMVERTICES2D];
198+ IssmDouble lsf[NUMVERTICES];
199+ int* indices = xNew<int>(NUMVERTICES2D);
200+
201+ /*Retrieve all inputs and parameters*/
202+ GetInputListOnVertices(&lsf[0],levelset_enum);
203+
204+ /* Determine distribution of ice over element.
205+ * Exploit: ice/no-ice parts are connected, so find starting vertex of segment*/
206+ lastindex=0;
207+ for(i=0;i<NUMVERTICES2D;i++){ // go backwards along vertices, and check for sign change
208+ ind0=(NUMVERTICES2D-i)%NUMVERTICES2D;
209+ ind1=(ind0-1+NUMVERTICES2D)%NUMVERTICES2D;
210+ if((lsf[ind0]-level)*(lsf[ind1]-level)<=0.){ // levelset has been crossed, find last index belonging to segment
211+ if(lsf[ind1]==level) //if levelset intersects 2nd vertex, choose this vertex as last
212+ lastindex=ind1;
213+ else
214+ lastindex=ind0;
215+ break;
216+ }
217+ }
218+
219+ numiceverts=0;
220+ numnoiceverts=0;
221+ for(i=0;i<NUMVERTICES2D;i++){
222+ ind0=(lastindex+i)%NUMVERTICES2D;
223+ if(lsf[i]<=level){
224+ indices_ice[numiceverts]=i;
225+ numiceverts++;
226+ }
227+ else{
228+ indices_noice[numnoiceverts]=i;
229+ numnoiceverts++;
230+ }
231+ }
232+ //merge indices
233+ for(i=0;i<numiceverts;i++){indices[i]=indices_ice[i];}
234+ for(i=0;i<numnoiceverts;i++){indices[numiceverts+i]=indices_noice[i];}
235+
236+ switch (numiceverts){
237+ case 0: // no vertex has ice: element is ice free, no intersection
238+ for(i=0;i<2;i++)
239+ fraction[i]=0.;
240+ break;
241+ case 1: // one vertex has ice:
242+ for(i=0;i<2;i++){
243+ fraction[i]=(level-lsf[indices[0]])/(lsf[indices[numiceverts+i]]-lsf[indices[0]]);
244+ }
245+ break;
246+ case 2: // two vertices have ice: fraction is computed from first ice vertex to last in CCW fashion
247+ for(i=0;i<2;i++){
248+ fraction[i]=(level-lsf[indices[i]])/(lsf[indices[numiceverts]]-lsf[indices[i]]);
249+ }
250+ break;
251+ case NUMVERTICES2D: // all vertices have ice: return triangle area
252+ for(i=0;i<2;i++)
253+ fraction[i]=1.;
254+ break;
255+ default:
256+ _error_("Wrong number of ice vertices in Penta::GetLevelsetIntersectionBase!");
257+ break;
258+ }
259+
260+ *pindices=indices;
261+ *pnumiceverts=numiceverts;
262+}
263+/*}}}*/
264 Node* Penta::GetNode(int node_number){/*{{{*/
265 _assert_(node_number>=0);
266 _assert_(node_number<this->NumberofNodes(this->element_type));
267@@ -2311,6 +2488,65 @@
268
269 }
270 /*}}}*/
271+void Penta::RignotMeltParameterization(){/*{{{*/
272+
273+ IssmDouble A, B, alpha, beta;
274+ IssmDouble bed,qsg,qsg_basin,TF,yts;
275+ int numbasins;
276+ IssmDouble basinid[NUMVERTICES];
277+ IssmDouble* basin_icefront_area=NULL;
278+
279+ /* Coefficients */
280+ A = 3e-4;
281+ B = 0.15;
282+ alpha = 0.39;
283+ beta = 1.18;
284+
285+ /*Get inputs*/
286+ Input* bed_input = this->GetInput(BedEnum); _assert_(bed_input);
287+ Input* qsg_input = this->GetInput(FrontalForcingsSubglacialDischargeEnum); _assert_(qsg_input);
288+ Input* TF_input = this->GetInput(FrontalForcingsThermalForcingEnum); _assert_(TF_input);
289+ GetInputListOnVertices(&basinid[0],FrontalForcingsBasinIdEnum);
290+
291+ this->FindParam(&yts, ConstantsYtsEnum);
292+ this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
293+ this->parameters->FindParam(&basin_icefront_area,&numbasins,FrontalForcingsBasinIcefrontAreaEnum);
294+
295+ IssmDouble meltrates[NUMVERTICES2D]; //frontal melt-rate
296+
297+ /* Start looping on the number of vertices: */
298+ GaussPenta* gauss=new GaussPenta();
299+ for(int iv=0;iv<NUMVERTICES2D;iv++){
300+ gauss->GaussVertex(iv);
301+
302+ /* Get variables */
303+ bed_input->GetInputValue(&bed,gauss);
304+ qsg_input->GetInputValue(&qsg,gauss);
305+ TF_input->GetInputValue(&TF,gauss);
306+
307+ if(basin_icefront_area[reCast<int>(basinid[iv])-1]==0.) meltrates[iv]=0.;
308+ else{
309+ /* change the unit of qsg (m^3/d -> m/d) with ice front area */
310+ qsg_basin=qsg/basin_icefront_area[reCast<int>(basinid[iv])-1];
311+
312+ /* calculate melt rates */
313+ meltrates[iv]=(A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta)/86400; //[m/s]
314+ }
315+
316+ if(xIsNan<IssmDouble>(meltrates[iv])) _error_("NaN found in vector");
317+ if(xIsInf<IssmDouble>(meltrates[iv])) _error_("Inf found in vector");
318+ }
319+
320+ /*Add input*/
321+ this->inputs->AddInput(new PentaInput(CalvingMeltingrateEnum,&meltrates[0],P1Enum));
322+
323+ this->InputExtrude(CalvingMeltingrateEnum,-1);
324+
325+ /*Cleanup and return*/
326+ xDelete<IssmDouble>(basin_icefront_area);
327+ delete gauss;
328+}
329+/*}}}*/
330 void Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int N, int M){/*{{{*/
331
332 IssmDouble values[NUMVERTICES];
333Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
334===================================================================
335--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23800)
336+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23801)
337@@ -3419,8 +3419,8 @@
338 IssmDouble* basin_icefront_area=NULL;
339
340 /* Coefficients */
341- A = 3e-4; //
342- B = 0.15; //
343+ A = 3e-4;
344+ B = 0.15;
345 alpha = 0.39;
346 beta = 1.18;
347
348@@ -3457,7 +3457,6 @@
349
350 if(xIsNan<IssmDouble>(meltrates[iv])) _error_("NaN found in vector");
351 if(xIsInf<IssmDouble>(meltrates[iv])) _error_("Inf found in vector");
352-
353 }
354
355 /*Add input*/
Note: See TracBrowser for help on using the repository browser.