source: issm/oecreview/Archive/20545-21336/ISSM-20718-20719.diff

Last change on this file was 21337, checked in by Mathieu Morlighem, 8 years ago

CHG: added Archive/20545-21336

File size: 7.0 KB
RevLine 
[21337]1Index: ../trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp (revision 20718)
4+++ ../trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp (revision 20719)
5@@ -26,6 +26,9 @@
6 femmodel->CalvingRateDevx();
7 femmodel->ElementOperationx(&Element::CalvingRateDev);
8 break;
9+ case CalvingMinthicknessEnum:
10+ femmodel->CalvingRateMinthicknessx();
11+ break;
12 default:
13 _error_("Caving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
14 }
15Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
16===================================================================
17--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 20718)
18+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 20719)
19@@ -2076,6 +2076,9 @@
20 case CalvingDevEnum:
21 this->CalvingRateDev();
22 break;
23+ case CalvingMinthicknessEnum:
24+ this->CalvingRateMinthickness();
25+ break;
26 default:
27 _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
28 }
29Index: ../trunk-jpl/src/c/classes/Elements/Element.h
30===================================================================
31--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 20718)
32+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 20719)
33@@ -175,6 +175,7 @@
34 virtual void AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0;
35 virtual void CalvingRateLevermann(void)=0;
36 virtual void CalvingRateDev(void){_error_("not implemented yet");};
37+ virtual void CalvingRateMinthickness(void){_error_("not implemented yet");};
38 virtual void WriteLevelsetSegment(DataSet* segments){_error_("not implemented yet");};
39 virtual void ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){_error_("not implemented yet");};
40 virtual IssmDouble CharacteristicLength(void)=0;
41Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
42===================================================================
43--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 20718)
44+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 20719)
45@@ -337,6 +337,51 @@
46 delete gauss;
47 }
48 /*}}}*/
49+void Tria::CalvingRateMinthickness(){/*{{{*/
50+
51+ IssmDouble calvingratex[NUMVERTICES];
52+ IssmDouble calvingratey[NUMVERTICES];
53+ IssmDouble calvingrate[NUMVERTICES];
54+ IssmDouble vx,vy,H,minH;
55+
56+ /*Retrieve all inputs and parameters we will need*/
57+ Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
58+ Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
59+ Input* H_input = inputs->GetInput(ThicknessEnum); _assert_(H_input);
60+ parameters->FindParam(&minH,CalvingMinthicknessEnum);
61+
62+ /* Start looping on the number of vertices: */
63+ GaussTria* gauss=new GaussTria();
64+ for(int iv=0;iv<NUMVERTICES;iv++){
65+ gauss->GaussVertex(iv);
66+
67+ /*Get velocity components and thickness*/
68+ vx_input->GetInputValue(&vx,gauss);
69+ vy_input->GetInputValue(&vy,gauss);
70+ H_input->GetInputValue(&H,gauss);
71+
72+ /*Assign values*/
73+ if(H>minH){
74+ calvingratex[iv]=0.;
75+ calvingratey[iv]=0.;
76+ calvingrate[iv]=0.;
77+ }
78+ else{
79+ calvingratex[iv]=vx+1e-2; //counter balance advance + add some retreat
80+ calvingratey[iv]=vy+1e-2;
81+ calvingrate[iv]=sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
82+ }
83+ }
84+
85+ /*Add input*/
86+ this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
87+ this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
88+ this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
89+
90+ /*Clean up and return*/
91+ delete gauss;
92+}
93+/*}}}*/
94 void Tria::WriteLevelsetSegment(DataSet* segments){/*{{{*/
95
96 if(!this->IsZeroLevelset(MaskIceLevelsetEnum)) return;
97Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
98===================================================================
99--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 20718)
100+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 20719)
101@@ -53,6 +53,7 @@
102 void AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
103 void CalvingRateLevermann();
104 void CalvingRateDev();
105+ void CalvingRateMinthickness();
106 void WriteLevelsetSegment(DataSet* segments);
107 void ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments);
108 IssmDouble CharacteristicLength(void);
109Index: ../trunk-jpl/src/c/classes/FemModel.h
110===================================================================
111--- ../trunk-jpl/src/c/classes/FemModel.h (revision 20718)
112+++ ../trunk-jpl/src/c/classes/FemModel.h (revision 20719)
113@@ -94,6 +94,7 @@
114 void DeviatoricStressx();
115 void CalvingRateLevermannx();
116 void CalvingRateDevx();
117+ void CalvingRateMinthicknessx();
118 void ResetLevelset();
119 #ifdef _HAVE_DAKOTA_
120 void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);
121Index: ../trunk-jpl/src/c/classes/FemModel.cpp
122===================================================================
123--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 20718)
124+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 20719)
125@@ -2165,6 +2165,14 @@
126 }
127 }
128 /*}}}*/
129+void FemModel::CalvingRateMinthicknessx(){/*{{{*/
130+
131+ for(int i=0;i<elements->Size();i++){
132+ Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
133+ element->CalvingRateMinthickness();
134+ }
135+}
136+/*}}}*/
137 void FemModel::StrainRateparallelx(){/*{{{*/
138
139 for(int i=0;i<elements->Size();i++){
140Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
141===================================================================
142--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 20718)
143+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 20719)
144@@ -197,6 +197,7 @@
145 switch(calvinglaw){
146 case DefaultCalvingEnum:
147 case CalvingDevEnum:
148+ case CalvingMinthicknessEnum:
149 lsf_slopex_input = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
150 if(dim==2) lsf_slopey_input = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
151 calvingrate_input = basalelement->GetInput(CalvingCalvingrateEnum); _assert_(calvingrate_input);
152@@ -250,6 +251,8 @@
153 /*Get calving speed*/
154 switch(calvinglaw){
155 case DefaultCalvingEnum:
156+ case CalvingDevEnum:
157+ case CalvingMinthicknessEnum:
158 lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
159 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
160 calvingrate_input->GetInputValue(&calvingrate,gauss);
161@@ -270,7 +273,6 @@
162 break;
163
164 case CalvingLevermannEnum:
165- case CalvingDevEnum:
166 calvingratex_input->GetInputValue(&c[0],gauss);
167 if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
168 meltingrate_input->GetInputValue(&meltingrate,gauss);
Note: See TracBrowser for help on using the repository browser.