source: issm/oecreview/Archive/18296-19100/ISSM-18757-18758.diff@ 19102

Last change on this file since 19102 was 19102, checked in by Mathieu Morlighem, 10 years ago

NEW: added 18296-19100
File size: 28.4 KB
RevLine 
[19102]1Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 18757)
4+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 18758)
5@@ -19,10 +19,9 @@
6 }
7 /*}}}*/
8 void LevelsetAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
9- int finiteelement;
10
11 /*Finite element type*/
12- finiteelement = P1Enum;
13+ int finiteelement = P1Enum;
14
15 /*Update elements: */
16 int counter=0;
17@@ -37,7 +36,20 @@
18 iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
19 iomodel->FetchDataToInput(elements,VxEnum);
20 iomodel->FetchDataToInput(elements,VyEnum);
21- iomodel->FetchDataToInput(elements,MasstransportCalvingrateEnum);
22+
23+ /*Get calving parameters*/
24+ int calvinglaw;
25+ iomodel->Constant(&calvinglaw,CalvingLawEnum);
26+ switch(calvinglaw){
27+ case DefaultCalvingEnum:
28+ iomodel->FetchDataToInput(elements,CalvingCalvingrateEnum);
29+ break;
30+ case CalvingLevermannEnum:
31+ iomodel->FetchDataToInput(elements,CalvinglevermannCoeffEnum);
32+ break;
33+ default:
34+ _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
35+ }
36 }
37 /*}}}*/
38 void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
39@@ -99,19 +111,21 @@
40 Element* basalelement = element->SpawnBasalElement();
41
42 /*Intermediaries */
43- int dim, domaintype;
44- bool iscalvingrate;
45+ int stabilization=2;
46+ int dim, domaintype, calvinglaw;
47+ bool iscalving;
48 int i, row, col;
49 IssmDouble kappa;
50 IssmDouble Jdet, dt, D_scalar;
51 IssmDouble h,hx,hy,hz;
52 IssmDouble vel;
53-// IssmDouble norm_dlsf;
54+ IssmDouble norm_dlsf, calvingrate;
55 IssmDouble* xyz_list = NULL;
56
57 /*Get problem dimension and whether there is calving or not*/
58- basalelement->FindParam(&iscalvingrate,MasstransportIscalvingrateEnum);
59+ basalelement->FindParam(&iscalving,TransientIscalvingEnum);
60 basalelement->FindParam(&domaintype,DomainTypeEnum);
61+ basalelement->FindParam(&calvinglaw,CalvingLawEnum);
62 switch(domaintype){
63 case Domain2DverticalEnum: dim = 1; break;
64 case Domain2DhorizontalEnum: dim = 2; break;
65@@ -130,45 +144,63 @@
66 IssmDouble* D = xNew<IssmDouble>(dim*dim);
67 IssmDouble* v = xNew<IssmDouble>(dim);
68 IssmDouble* w = xNew<IssmDouble>(dim);
69- IssmDouble* c = xNew<IssmDouble>(dim);
70- //IssmDouble* dlsf = xNew<IssmDouble>(dim);
71+ IssmDouble* c = xNewZeroInit<IssmDouble>(dim);
72+ IssmDouble* dlsf = xNew<IssmDouble>(dim);
73
74 /*Retrieve all inputs and parameters*/
75 basalelement->GetVerticesCoordinates(&xyz_list);
76 basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
77- Input* vx_input=NULL;
78- Input* vy_input=NULL;
79- Input* calvingratex_input=NULL;
80- Input* calvingratey_input=NULL;
81+ Input* vx_input = NULL;
82+ Input* vy_input = NULL;
83+ Input* calvingratex_input = NULL;
84+ Input* calvingratey_input = NULL;
85+ Input* lsf_slopex_input = NULL;
86+ Input* lsf_slopey_input = NULL;
87+ Input* calvingrate_input = NULL;
88+
89+ /*Load velocities*/
90 if(domaintype==Domain2DhorizontalEnum){
91 vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
92 vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
93- if(iscalvingrate){
94- calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
95- calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
96- }
97 }
98 else{
99 if(dim==1){
100 vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
101- if(iscalvingrate){
102- calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
103- }
104 }
105 if(dim==2){
106 vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
107 vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
108- if(iscalvingrate){
109- calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
110- calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
111- }
112 }
113 }
114
115-// Input* lsf_slopex_input = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
116-// Input* lsf_slopey_input = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
117- //Input* calvingrate_input = basalelement->GetInput(MasstransportCalvingrateEnum); _assert_(calvingrate_input);
118-
119+ /*Load calving inputs*/
120+ if(iscalving){
121+ switch(calvinglaw){
122+ case DefaultCalvingEnum:
123+ lsf_slopex_input = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
124+ if(dim==2) lsf_slopey_input = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
125+ calvingrate_input = basalelement->GetInput(CalvingCalvingrateEnum); _assert_(calvingrate_input);
126+ break;
127+ case CalvingLevermannEnum:
128+ if(domaintype==Domain2DhorizontalEnum){
129+ calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
130+ calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
131+ }
132+ else{
133+ if(dim==1){
134+ calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
135+ }
136+ if(dim==2){
137+ calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
138+ calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
139+ }
140+ }
141+ break;
142+ default:
143+ _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
144+ }
145+ }
146+
147 /* Start looping on the number of gaussian points: */
148 Gauss* gauss=basalelement->NewGauss(2);
149 for(int ig=gauss->begin();ig<gauss->end();ig++){
150@@ -192,34 +224,45 @@
151 GetBprime(Bprime,basalelement,xyz_list,gauss);
152 vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity
153 vy_input->GetInputValue(&v[1],gauss);
154- if(iscalvingrate){
155- calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity
156- calvingratey_input->GetInputValue(&c[1],gauss);
157- for(i=0;i<dim;i++) w[i]=v[i]-c[i];
158+
159+ /*Get calving speed*/
160+ if(iscalving){
161+ switch(calvinglaw){
162+ case DefaultCalvingEnum:
163+ lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
164+ if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
165+ calvingrate_input->GetInputValue(&calvingrate,gauss);
166+
167+ norm_dlsf=0.;
168+ for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
169+ norm_dlsf=sqrt(norm_dlsf);
170+
171+ if(norm_dlsf>1.e-10)
172+ for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
173+ else
174+ for(i=0;i<dim;i++) c[i]=0.;
175+ break;
176+ case CalvingLevermannEnum:
177+ calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity
178+ if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
179+ break;
180+ default:
181+ _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
182+ }
183 }
184- else{
185- for(i=0;i<dim;i++) w[i]=v[i];
186- }
187- //lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
188- //lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
189- //calvingrate_input->GetInputValue(&calvingrate,gauss);
190
191- //norm_dlsf=0.;
192- //for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
193- //norm_dlsf=sqrt(norm_dlsf);
194+ /*Levelset speed is ice velocity - calving rate*/
195+ for(i=0;i<dim;i++) w[i]=v[i]-c[i];
196
197- //if(norm_dlsf>1.e-10)
198- // for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
199- //else
200- //for(i=0;i<dim;i++) c[i]=0.;
201-
202-
203- for(row=0;row<dim;row++)
204- for(col=0;col<dim;col++)
205+ /*Compute D*/
206+ for(row=0;row<dim;row++){
207+ for(col=0;col<dim;col++){
208 if(row==col)
209- D[row*dim+col]=D_scalar*w[row];
210+ D[row*dim+col]=D_scalar*w[row];
211 else
212- D[row*dim+col]=0.;
213+ D[row*dim+col]=0.;
214+ }
215+ }
216
217 TripleMultiply(B,dim,numnodes,1,
218 D,dim,dim,0,
219@@ -227,9 +270,8 @@
220 &Ke->values[0],1);
221
222 /* Stabilization */
223- int stabilization=2;
224 vel=0.;
225- for(i=0;i<dim;i++) vel+=w[i]*w[i];
226+ for(i=0;i<dim;i++) vel+=v[i]*v[i];
227 vel=sqrt(vel)+1.e-14;
228 switch(stabilization){
229 case 0:
230@@ -255,10 +297,10 @@
231 case 2:
232 /* Streamline Upwinding */
233 basalelement->ElementSizes(&hx,&hy,&hz);
234- h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
235+ h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) );
236 for(row=0;row<dim;row++)
237 for(col=0;col<dim;col++)
238- D[row*dim+col] = D_scalar*h/(2.*vel)*w[row]*w[col];
239+ D[row*dim+col] = D_scalar*h/(2.*vel)*v[row]*v[col];
240
241 TripleMultiply(Bprime,dim,numnodes,1,
242 D,dim,dim,0,
243@@ -279,7 +321,7 @@
244 xDelete<IssmDouble>(v);
245 xDelete<IssmDouble>(w);
246 xDelete<IssmDouble>(c);
247- //xDelete<IssmDouble>(dlsf);
248+ xDelete<IssmDouble>(dlsf);
249 delete gauss;
250 if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
251 return Ke;
252Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
253===================================================================
254--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18757)
255+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18758)
256@@ -208,6 +208,12 @@
257 DamageEnum,
258 NewDamageEnum,
259 StressIntensityFactorEnum,
260+ CalvingLawEnum,
261+ CalvingCalvingrateEnum,
262+ CalvingLevermannEnum,
263+ DefaultCalvingEnum,
264+ CalvingRequestedOutputsEnum,
265+ CalvinglevermannCoeffEnum,
266 CalvingratexEnum,
267 CalvingrateyEnum,
268 CalvingratexAverageEnum,
269@@ -249,13 +255,10 @@
270 Domain3DEnum,
271 MiscellaneousNameEnum, //FIXME: only used by qmu, should not be marshalled (already in queueing script)
272 MasstransportHydrostaticAdjustmentEnum,
273- MasstransportIscalvingrateEnum,
274- MasstransportLevermannCalvingCoeffEnum,
275 MasstransportIsfreesurfaceEnum,
276 MasstransportMinThicknessEnum,
277 MasstransportPenaltyFactorEnum,
278 MasstransportSpcthicknessEnum,
279- MasstransportCalvingrateEnum,
280 MasstransportStabilizationEnum,
281 MasstransportVertexPairingEnum,
282 MasstransportNumRequestedOutputsEnum,
283@@ -313,6 +316,7 @@
284 TransientIsgiaEnum,
285 TransientIsdamageevolutionEnum,
286 TransientIshydrologyEnum,
287+ TransientIscalvingEnum,
288 TransientNumRequestedOutputsEnum,
289 TransientRequestedOutputsEnum,
290 PotentialEnum,
291Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
292===================================================================
293--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18757)
294+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18758)
295@@ -216,6 +216,12 @@
296 case DamageEnum : return "Damage";
297 case NewDamageEnum : return "NewDamage";
298 case StressIntensityFactorEnum : return "StressIntensityFactor";
299+ case CalvingLawEnum : return "CalvingLaw";
300+ case CalvingCalvingrateEnum : return "CalvingCalvingrate";
301+ case CalvingLevermannEnum : return "CalvingLevermann";
302+ case DefaultCalvingEnum : return "DefaultCalving";
303+ case CalvingRequestedOutputsEnum : return "CalvingRequestedOutputs";
304+ case CalvinglevermannCoeffEnum : return "CalvinglevermannCoeff";
305 case CalvingratexEnum : return "Calvingratex";
306 case CalvingrateyEnum : return "Calvingratey";
307 case CalvingratexAverageEnum : return "CalvingratexAverage";
308@@ -257,13 +263,10 @@
309 case Domain3DEnum : return "Domain3D";
310 case MiscellaneousNameEnum : return "MiscellaneousName";
311 case MasstransportHydrostaticAdjustmentEnum : return "MasstransportHydrostaticAdjustment";
312- case MasstransportIscalvingrateEnum : return "MasstransportIscalvingrate";
313- case MasstransportLevermannCalvingCoeffEnum : return "MasstransportLevermannCalvingCoeff";
314 case MasstransportIsfreesurfaceEnum : return "MasstransportIsfreesurface";
315 case MasstransportMinThicknessEnum : return "MasstransportMinThickness";
316 case MasstransportPenaltyFactorEnum : return "MasstransportPenaltyFactor";
317 case MasstransportSpcthicknessEnum : return "MasstransportSpcthickness";
318- case MasstransportCalvingrateEnum : return "MasstransportCalvingrate";
319 case MasstransportStabilizationEnum : return "MasstransportStabilization";
320 case MasstransportVertexPairingEnum : return "MasstransportVertexPairing";
321 case MasstransportNumRequestedOutputsEnum : return "MasstransportNumRequestedOutputs";
322@@ -321,6 +324,7 @@
323 case TransientIsgiaEnum : return "TransientIsgia";
324 case TransientIsdamageevolutionEnum : return "TransientIsdamageevolution";
325 case TransientIshydrologyEnum : return "TransientIshydrology";
326+ case TransientIscalvingEnum : return "TransientIscalving";
327 case TransientNumRequestedOutputsEnum : return "TransientNumRequestedOutputs";
328 case TransientRequestedOutputsEnum : return "TransientRequestedOutputs";
329 case PotentialEnum : return "Potential";
330Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
331===================================================================
332--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18757)
333+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18758)
334@@ -219,6 +219,12 @@
335 else if (strcmp(name,"Damage")==0) return DamageEnum;
336 else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
337 else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
338+ else if (strcmp(name,"CalvingLaw")==0) return CalvingLawEnum;
339+ else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
340+ else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
341+ else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
342+ else if (strcmp(name,"CalvingRequestedOutputs")==0) return CalvingRequestedOutputsEnum;
343+ else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
344 else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
345 else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
346 else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
347@@ -253,23 +259,20 @@
348 else if (strcmp(name,"MeshY")==0) return MeshYEnum;
349 else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
350 else if (strcmp(name,"MeshElementtype")==0) return MeshElementtypeEnum;
351- else if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
352+ else stage=3;
353+ }
354+ if(stage==3){
355+ if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
356 else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum;
357 else if (strcmp(name,"Domain2Dhorizontal")==0) return Domain2DhorizontalEnum;
358 else if (strcmp(name,"Domain2Dvertical")==0) return Domain2DverticalEnum;
359 else if (strcmp(name,"Domain3D")==0) return Domain3DEnum;
360 else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum;
361- else stage=3;
362- }
363- if(stage==3){
364- if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum;
365- else if (strcmp(name,"MasstransportIscalvingrate")==0) return MasstransportIscalvingrateEnum;
366- else if (strcmp(name,"MasstransportLevermannCalvingCoeff")==0) return MasstransportLevermannCalvingCoeffEnum;
367+ else if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum;
368 else if (strcmp(name,"MasstransportIsfreesurface")==0) return MasstransportIsfreesurfaceEnum;
369 else if (strcmp(name,"MasstransportMinThickness")==0) return MasstransportMinThicknessEnum;
370 else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum;
371 else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum;
372- else if (strcmp(name,"MasstransportCalvingrate")==0) return MasstransportCalvingrateEnum;
373 else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
374 else if (strcmp(name,"MasstransportVertexPairing")==0) return MasstransportVertexPairingEnum;
375 else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum;
376@@ -327,6 +330,7 @@
377 else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
378 else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum;
379 else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum;
380+ else if (strcmp(name,"TransientIscalving")==0) return TransientIscalvingEnum;
381 else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum;
382 else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
383 else if (strcmp(name,"Potential")==0) return PotentialEnum;
384@@ -378,14 +382,14 @@
385 else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
386 else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
387 else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
388- else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
389+ else stage=4;
390+ }
391+ if(stage==4){
392+ if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
393 else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
394 else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
395 else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
396- else stage=4;
397- }
398- if(stage==4){
399- if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
400+ else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
401 else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
402 else if (strcmp(name,"BalancethicknessSoftAnalysis")==0) return BalancethicknessSoftAnalysisEnum;
403 else if (strcmp(name,"BalancethicknessSoftSolution")==0) return BalancethicknessSoftSolutionEnum;
404@@ -501,14 +505,14 @@
405 else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
406 else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
407 else if (strcmp(name,"StringParam")==0) return StringParamEnum;
408- else if (strcmp(name,"Seg")==0) return SegEnum;
409+ else stage=5;
410+ }
411+ if(stage==5){
412+ if (strcmp(name,"Seg")==0) return SegEnum;
413 else if (strcmp(name,"SegInput")==0) return SegInputEnum;
414 else if (strcmp(name,"Tria")==0) return TriaEnum;
415 else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
416- else stage=5;
417- }
418- if(stage==5){
419- if (strcmp(name,"Tetra")==0) return TetraEnum;
420+ else if (strcmp(name,"Tetra")==0) return TetraEnum;
421 else if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
422 else if (strcmp(name,"Penta")==0) return PentaEnum;
423 else if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
424@@ -624,14 +628,14 @@
425 else if (strcmp(name,"P2")==0) return P2Enum;
426 else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
427 else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum;
428- else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
429+ else stage=6;
430+ }
431+ if(stage==6){
432+ if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
433 else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
434 else if (strcmp(name,"P1xP3")==0) return P1xP3Enum;
435 else if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
436- else stage=6;
437- }
438- if(stage==6){
439- if (strcmp(name,"P1P1")==0) return P1P1Enum;
440+ else if (strcmp(name,"P1P1")==0) return P1P1Enum;
441 else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
442 else if (strcmp(name,"MINI")==0) return MINIEnum;
443 else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
444@@ -747,14 +751,14 @@
445 else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
446 else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
447 else if (strcmp(name,"TransientIslevelset")==0) return TransientIslevelsetEnum;
448- else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
449+ else stage=7;
450+ }
451+ if(stage==7){
452+ if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
453 else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
454 else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
455 else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
456- else stage=7;
457- }
458- if(stage==7){
459- if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
460+ else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
461 else if (strcmp(name,"Seaiceatm")==0) return SeaiceatmEnum;
462 else if (strcmp(name,"Seaiceocean")==0) return SeaiceoceanEnum;
463 else if (strcmp(name,"SeaiceThickness")==0) return SeaiceThicknessEnum;
464Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
465===================================================================
466--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 18757)
467+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 18758)
468@@ -66,6 +66,7 @@
469 parameters->AddObject(iomodel->CopyConstantObject(QmuIsdakotaEnum));
470 parameters->AddObject(iomodel->CopyConstantObject(InversionIscontrolEnum));
471 parameters->AddObject(iomodel->CopyConstantObject(InversionTypeEnum));
472+ parameters->AddObject(iomodel->CopyConstantObject(CalvingLawEnum));
473 if(solution_type==SeaiceSolutionEnum){
474 }
475 else{
476@@ -82,8 +83,7 @@
477 parameters->AddObject(iomodel->CopyConstantObject(TransientIslevelsetEnum));
478 parameters->AddObject(iomodel->CopyConstantObject(TransientIsdamageevolutionEnum));
479 parameters->AddObject(iomodel->CopyConstantObject(TransientIshydrologyEnum));
480- parameters->AddObject(iomodel->CopyConstantObject(MasstransportIscalvingrateEnum));
481- parameters->AddObject(iomodel->CopyConstantObject(MasstransportLevermannCalvingCoeffEnum));
482+ parameters->AddObject(iomodel->CopyConstantObject(TransientIscalvingEnum));
483 parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
484 parameters->AddObject(iomodel->CopyConstantObject(GiaCrossSectionShapeEnum));
485
486Index: ../trunk-jpl/src/c/cores/transient_core.cpp
487===================================================================
488--- ../trunk-jpl/src/c/cores/transient_core.cpp (revision 18757)
489+++ ../trunk-jpl/src/c/cores/transient_core.cpp (revision 18758)
490@@ -20,11 +20,11 @@
491
492 /*parameters: */
493 IssmDouble starttime,finaltime,dt,yts;
494- bool isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology,iscalvingrate;
495+ bool isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology,iscalving;
496 bool save_results,dakota_analysis;
497 bool time_adapt=false;
498 int output_frequency;
499- int domaintype,groundingline_migration;
500+ int domaintype,groundingline_migration,calvinglaw;
501 int numoutputs = 0;
502 Analysis *analysis = NULL;
503 char** requested_outputs = NULL;
504@@ -52,7 +52,8 @@
505 femmodel->parameters->FindParam(&isdamageevolution,TransientIsdamageevolutionEnum);
506 femmodel->parameters->FindParam(&ishydrology,TransientIshydrologyEnum);
507 femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
508- femmodel->parameters->FindParam(&iscalvingrate,MasstransportIscalvingrateEnum);
509+ femmodel->parameters->FindParam(&iscalving,TransientIscalvingEnum);
510+ femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
511 if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
512 femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
513 if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum);
514@@ -101,7 +102,7 @@
515 }
516
517 if(islevelset){
518- if(iscalvingrate){
519+ if(iscalving && calvinglaw==CalvingLevermannEnum){
520 if(VerboseSolution()) _printf0_(" computing calving rate\n");
521 femmodel->StrainRateparallelx();
522 femmodel->StrainRateperpendicularx();
523Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
524===================================================================
525--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18757)
526+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18758)
527@@ -438,7 +438,7 @@
528 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
529 Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum); _assert_(strainparallel_input);
530 Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum); _assert_(strainperpendicular_input);
531- this->parameters->FindParam(&propcoeff,MasstransportLevermannCalvingCoeffEnum);
532+ Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum); _assert_(levermanncoeff_input);
533
534 /* Start looping on the number of vertices: */
535 gauss=new GaussTria();
536@@ -451,6 +451,7 @@
537 vel=vx*vx+vy*vy;
538 strainparallel_input->GetInputValue(&strainparallel,gauss);
539 strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
540+ levermanncoeff_input->GetInputValue(&propcoeff,gauss);
541
542 /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
543 calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
544@@ -464,7 +465,7 @@
545 /*Add input*/
546 this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
547 this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
548- this->inputs->AddInput(new TriaInput(MasstransportCalvingrateEnum,&calvingrate[0],P1Enum));
549+ this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
550
551 /*Clean up and return*/
552 delete gauss;
553Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
554===================================================================
555--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18757)
556+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18758)
557@@ -366,8 +366,8 @@
558 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
559 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
560 Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum); _assert_(strainparallel_input);
561- Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum); _assert_(strainperpendicular_input);
562- this->parameters->FindParam(&propcoeff,MasstransportLevermannCalvingCoeffEnum);
563+ Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum); _assert_(strainperpendicular_input);
564+ Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum); _assert_(levermanncoeff_input);
565
566 /* Start looping on the number of vertices: */
567 gauss=new GaussPenta();
568@@ -380,6 +380,7 @@
569 vel=vx*vx+vy*vy;
570 strainparallel_input->GetInputValue(&strainparallel,gauss);
571 strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
572+ levermanncoeff_input->GetInputValue(&propcoeff,gauss);
573
574 /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
575 calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
576@@ -393,7 +394,7 @@
577 /*Add input*/
578 this->inputs->AddInput(new PentaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
579 this->inputs->AddInput(new PentaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
580- this->inputs->AddInput(new PentaInput(MasstransportCalvingrateEnum,&calvingrate[0],P1Enum));
581+ this->inputs->AddInput(new PentaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
582
583 /*Clean up and return*/
584 delete gauss;
585Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
586===================================================================
587--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18757)
588+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18758)
589@@ -1091,7 +1091,7 @@
590 name==LevelsetfunctionSlopeXEnum ||
591 name==LevelsetfunctionSlopeYEnum ||
592 name==LevelsetfunctionPicardEnum ||
593- name==MasstransportCalvingrateEnum ||
594+ //name==CalvingCalvingrateEnum ||
595 name==GradientEnum ||
596 name==OldGradientEnum ||
597 name==ConvergedEnum ||
598@@ -1211,7 +1211,7 @@
599 break;
600 case CalvingratexEnum:
601 case CalvingrateyEnum:
602- case MasstransportCalvingrateEnum:
603+ case CalvingCalvingrateEnum:
604 this->StrainRateparallel();
605 this->StrainRateperpendicular();
606 this->CalvingRateLevermann();
Note: See TracBrowser for help on using the repository browser.