source: issm/oecreview/Archive/23390-24306/ISSM-23606-23607.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: 10.2 KB
RevLine 
[24307]1Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
2===================================================================
3--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23606)
4+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23607)
5@@ -331,17 +331,17 @@
6 IssmDouble xyz_list[NUMVERTICES][3];
7 IssmDouble calvingrate[NUMVERTICES];
8 IssmDouble vx,vy,vel;
9- IssmDouble critical_fraction,water_height;
10- IssmDouble bed,Ho,thickness,float_depth;
11+ IssmDouble water_height, bed,Ho,thickness,surface;
12 IssmDouble surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal;
13- IssmDouble strainparallel, straineffective;
14+ IssmDouble B, strainparallel, straineffective;
15 IssmDouble s_xx,s_xy,s_yy,s1,s2,stmp;
16-
17+ int crevasse_opening_stress;
18+
19 /* Get node coordinates and dof list: */
20 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
21
22- /*Get the critical fraction of thickness surface and basal crevasses penetrate for calving onset*/
23- this->parameters->FindParam(&critical_fraction,CalvingCrevasseDepthEnum);
24+ /*retrieve the type of crevasse_opening_stress*/
25+ this->parameters->FindParam(&crevasse_opening_stress,CalvingCrevasseDepthEnum);
26
27 IssmDouble rho_ice = this->GetMaterialParameter(MaterialsRhoIceEnum);
28 IssmDouble rho_seawater = this->GetMaterialParameter(MaterialsRhoSeawaterEnum);
29@@ -360,6 +360,7 @@
30 Input* s_xx_input = inputs->GetInput(DeviatoricStressxxEnum); _assert_(s_xx_input);
31 Input* s_xy_input = inputs->GetInput(DeviatoricStressxyEnum); _assert_(s_xy_input);
32 Input* s_yy_input = inputs->GetInput(DeviatoricStressyyEnum); _assert_(s_yy_input);
33+ Input* B_input = inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input);
34
35 /*Loop over all elements of this partition*/
36 GaussTria* gauss=new GaussTria();
37@@ -368,7 +369,7 @@
38
39 H_input->GetInputValue(&thickness,gauss);
40 bed_input->GetInputValue(&bed,gauss);
41- surface_input->GetInputValue(&float_depth,gauss);
42+ surface_input->GetInputValue(&surface,gauss);
43 strainrateparallel_input->GetInputValue(&strainparallel,gauss);
44 strainrateeffective_input->GetInputValue(&straineffective,gauss);
45 vx_input->GetInputValue(&vx,gauss);
46@@ -377,6 +378,7 @@
47 s_xx_input->GetInputValue(&s_xx,gauss);
48 s_xy_input->GetInputValue(&s_xy,gauss);
49 s_yy_input->GetInputValue(&s_yy,gauss);
50+ B_input->GetInputValue(&B,gauss);
51
52 vel=sqrt(vx*vx+vy*vy)+1.e-14;
53
54@@ -387,28 +389,31 @@
55 Ho = thickness - (rho_seawater/rho_ice) * (-bed);
56 if(Ho<0.) Ho=0.;
57
58- /*Otero2010: balance between the tensile deviatoric stress and ice overburden pressure*/
59- /*surface crevasse*/
60- //surface_crevasse[iv] = rheology_B * strainparallel * pow(straineffective, ((1 / rheology_n)-1)) / (rho_ice * constant_g);
61- surface_crevasse[iv] = s1 / (rho_ice*constant_g);
62+ if(crevasse_opening_stress==0){ /*Otero2010: balance between the tensile deviatoric stress and ice overburden pressure*/
63+ surface_crevasse[iv] = B * strainparallel * pow(straineffective, ((1 / rheology_n)-1)) / (rho_ice * constant_g);
64+ basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (B * strainparallel * pow(straineffective,((1/rheology_n)-1)) / (rho_ice*constant_g) - Ho);
65+ }
66+ else if(crevasse_opening_stress==1){ /*Benn2017,Todd2018: maximum principal stress */
67+ surface_crevasse[iv] = s1 / (rho_ice*constant_g);
68+ basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho);
69+ }
70+
71+ /* some constraints */
72 if (surface_crevasse[iv]<0.) {
73 surface_crevasse[iv]=0.;
74 water_height = 0.;
75 }
76+ if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
77+ if (bed>0.) basal_crevasse[iv] = 0.;
78+
79 //if (surface_crevasse[iv]<water_height){
80 // water_height = surface_crevasse[iv];
81 //}
82-
83- /*basal crevasse*/
84- //basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (rheology_B * strainparallel * pow(straineffective,((1/rheology_n)-1)) / (rho_ice*constant_g) - Ho);
85- basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho);
86- if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
87- if (bed>0.) basal_crevasse[iv] = 0.;
88-
89- H_surf = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height - critical_fraction*float_depth;
90- H_surfbasal = (surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv])-(critical_fraction*thickness);
91-
92- crevasse_depth[iv] = max(H_surf,H_surfbasal);
93+
94+ /* add water in surface crevasse */
95+ surface_crevasse[iv] = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height; /* surface crevasse + water */
96+ crevasse_depth[iv] = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv]; /* surface crevasse + basal crevasse + water */
97+
98 }
99
100 this->inputs->AddInput(new TriaInput(SurfaceCrevasseEnum,&surface_crevasse[0],P1Enum));
101Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
102===================================================================
103--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 23606)
104+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 23607)
105@@ -106,7 +106,7 @@
106 case CalvingHabEnum:
107 break;
108 case CalvingCrevasseDepthEnum:
109- parameters->AddObject(iomodel->CopyConstantObject("md.calving.critical_fraction",CalvingCrevasseDepthEnum));
110+ parameters->AddObject(iomodel->CopyConstantObject("md.calving.crevasse_opening_stress",CalvingCrevasseDepthEnum));
111 break;
112 case CalvingDev2Enum:
113 parameters->AddObject(iomodel->CopyConstantObject("md.calving.height_above_floatation",CalvingHeightAboveFloatationEnum));
114@@ -670,7 +670,8 @@
115
116 /*Intermediaries*/
117 int calvinglaw;
118- IssmDouble min_thickness,thickness,hab_fraction,crevassedepth;
119+ IssmDouble min_thickness,thickness,hab_fraction;
120+ IssmDouble crevassedepth,surface_crevasse,surface,critical_fraction;
121 IssmDouble rho_ice,rho_water;
122 IssmDouble bed,water_depth;
123 IssmDouble levelset;
124@@ -757,16 +758,19 @@
125 /*Get the DistanceToCalvingfront*/
126 femmodel->elements->InputDuplicate(MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
127 femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
128-
129+
130 /*Vector of size number of nodes*/
131 vec_constraint_nodes=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes());
132
133 for(int i=0;i<femmodel->elements->Size();i++){
134- Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
135- int numnodes = element->GetNumberOfNodes();
136- Gauss* gauss = element->NewGauss();
137- Input* crevassedepth_input = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input);
138- Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input);
139+ Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
140+ int numnodes = element->GetNumberOfNodes();
141+ Gauss* gauss = element->NewGauss();
142+ Input* crevassedepth_input = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input);
143+ Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input);
144+ Input* surface_crevasse_input = element->GetInput(SurfaceCrevasseEnum); _assert_(surface_crevasse_input);
145+ Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
146+ Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
147
148 /*First, look at ice front and figure out if any of the nodes will be calved*/
149 if(element->IsIcefront()){
150@@ -775,7 +779,11 @@
151 Node* node=element->GetNode(in);
152 crevassedepth_input->GetInputValue(&crevassedepth,gauss);
153 bed_input->GetInputValue(&bed,gauss);
154- if(crevassedepth>0. && bed<0.){
155+ surface_crevasse_input->GetInputValue(&surface_crevasse,gauss);
156+ thickness_input->GetInputValue(&thickness,gauss);
157+ surface_input->GetInputValue(&surface,gauss);
158+
159+ if((surface_crevasse-surface>0. || crevassedepth-thickness>0.) && bed<0.){
160 vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL);
161 }
162 }
163@@ -790,12 +798,15 @@
164 while(nflipped){
165 local_nflipped=0;
166 for(int i=0;i<femmodel->elements->Size();i++){
167- Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
168- int numnodes = element->GetNumberOfNodes();
169- Gauss* gauss = element->NewGauss();
170- Input* levelset_input = element->GetInput(DistanceToCalvingfrontEnum); _assert_(levelset_input);
171- Input* crevassedepth_input = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input);
172- Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input);
173+ Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
174+ int numnodes = element->GetNumberOfNodes();
175+ Gauss* gauss = element->NewGauss();
176+ Input* levelset_input = element->GetInput(DistanceToCalvingfrontEnum); _assert_(levelset_input);
177+ Input* crevassedepth_input = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input);
178+ Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input);
179+ Input* surface_crevasse_input = element->GetInput(SurfaceCrevasseEnum); _assert_(surface_crevasse_input);
180+ Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
181+ Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
182
183 /*Is this element connected to a node that should be calved*/
184 bool isconnected = false;
185@@ -815,8 +826,11 @@
186 levelset_input->GetInputValue(&levelset,gauss);
187 crevassedepth_input->GetInputValue(&crevassedepth,gauss);
188 bed_input->GetInputValue(&bed,gauss);
189+ surface_crevasse_input->GetInputValue(&surface_crevasse,gauss);
190+ thickness_input->GetInputValue(&thickness,gauss);
191+ surface_input->GetInputValue(&surface,gauss);
192
193- if(crevassedepth>0. && bed<0. && levelset>-300. && levelset<0. && constraint_nodes[node->Sid()]==0.){
194+ if((surface_crevasse-surface>0. || crevassedepth-thickness>0.) && bed<0. && levelset>-300. && levelset<0. && constraint_nodes[node->Sid()]==0.){
195 local_nflipped++;
196 vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL);
197 }
Note: See TracBrowser for help on using the repository browser.