source: issm/oecreview/Archive/23390-24306/ISSM-23964-23965.diff

Last change on this file was 24307, checked in by Mathieu Morlighem, 5 years ago

NEW: adding Archive/23390-24306
File size: 8.5 KB
RevLine 
[24307]1Index: ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23964)
4+++ ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23965)
5@@ -371,7 +371,6 @@
6 }/*}}}*/
7 void HydrologyGlaDSAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
8 element->InputUpdateFromSolutionOneDof(solution,HydraulicPotentialEnum);
9- this->UpdateEffectivePressure(element);
10 }/*}}}*/
11 void HydrologyGlaDSAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
12 /*Default, do nothing*/
13@@ -509,3 +508,17 @@
14 delete gauss;
15 xDelete<IssmDouble>(N);
16 }/*}}}*/
17+void HydrologyGlaDSAnalysis::UpdateChannelCrossSection(FemModel* femmodel){/*{{{*/
18+
19+ bool ischannels;
20+ femmodel->parameters->FindParam(&ischannels,HydrologyIschannelsEnum);
21+ if(!ischannels) return;
22+
23+ for(int i=0;i<femmodel->loads->Size();i++){
24+ if(femmodel->loads->GetEnum(i)==ChannelEnum){
25+ Channel* channel=(Channel*)femmodel->loads->GetObjectByOffset(i);
26+ channel->UpdateChannelCrossSection();
27+ }
28+ }
29+
30+}/*}}}*/
31Index: ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h
32===================================================================
33--- ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h (revision 23964)
34+++ ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h (revision 23965)
35@@ -33,7 +33,8 @@
36 /*Specific to GlaDS*/
37 void UpdateSheetThickness(FemModel* femmodel);
38 void UpdateSheetThickness(Element* element);
39+ void UpdateChannelCrossSection(FemModel* femmodel);
40 void UpdateEffectivePressure(FemModel* femmodel);
41- void UpdateEffectivePressure(Element* element);
42+ void UpdateEffectivePressure(Element* element);
43 };
44 #endif
45Index: ../trunk-jpl/src/c/cores/hydrology_core.cpp
46===================================================================
47--- ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 23964)
48+++ ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 23965)
49@@ -197,8 +197,13 @@
50 InputDuplicatex(femmodel,HydraulicPotentialEnum,HydraulicPotentialOldEnum);
51 analysis->UpdateEffectivePressure(femmodel);
52 solutionsequence_shakti_nonlinear(femmodel);
53- if(VerboseSolution()) _printf0_(" updating sheet thickness\n");
54+
55+ if(VerboseSolution()) _printf0_(" updating effective pressure\n");
56+ analysis->UpdateEffectivePressure(femmodel);
57+ if(VerboseSolution()) _printf0_(" updating sheet thickness\n"); /*Uses N, so needs to come after*/
58 analysis->UpdateSheetThickness(femmodel);
59+ if(VerboseSolution()) _printf0_(" updating channels cross section\n");
60+ analysis->UpdateChannelCrossSection(femmodel);
61 delete analysis;
62 }
63
64Index: ../trunk-jpl/src/c/classes/Loads/Channel.cpp
65===================================================================
66--- ../trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 23964)
67+++ ../trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 23965)
68@@ -571,3 +571,126 @@
69 return pe;
70 }
71 /*}}}*/
72+void Channel::UpdateChannelCrossSection(void){/*{{{*/
73+
74+ /*Initialize Element matrix and return if necessary*/
75+ Tria* tria=(Tria*)element;
76+ if(!tria->IsIceInElement()){
77+ this->S = 0.;
78+ return;
79+ }
80+ _assert_(tria->FiniteElement()==P1Enum);
81+ int index1=tria->GetNodeIndex(nodes[0]);
82+ int index2=tria->GetNodeIndex(nodes[1]);
83+
84+ /*Intermediaries */
85+ IssmDouble Jdet,v1,qc,fFactor,Afactor,Bfactor,Xifactor;
86+ IssmDouble A,B,n,phi_old,phi,phi_0,dPw;
87+ IssmDouble H,h,b,dphi[2],dphids,dphimds,db[2],dbds;
88+ IssmDouble xyz_list[NUMVERTICES][3];
89+ IssmDouble xyz_list_tria[3][3];
90+ const int numnodes = NUMNODES;
91+
92+ /*Initialize Element vector and other vectors*/
93+ ElementMatrix* Ke=new ElementMatrix(this->nodes,NUMNODES,this->parameters);
94+ IssmDouble basis[NUMNODES];
95+ IssmDouble dbasisdx[2*NUMNODES];
96+ IssmDouble dbasisds[NUMNODES];
97+
98+ /*Retrieve all inputs and parameters*/
99+ GetVerticesCoordinates(&xyz_list[0][0] ,this->vertices,NUMVERTICES);
100+ GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3);
101+
102+ IssmDouble L = element->FindParam(MaterialsLatentheatEnum);
103+ IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
104+ IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
105+ IssmDouble g = element->FindParam(ConstantsGEnum);
106+ IssmDouble kc = element->FindParam(HydrologyChannelConductivityEnum);
107+ IssmDouble ks = element->FindParam(HydrologySheetConductivityEnum);
108+ IssmDouble lc = element->FindParam(HydrologyChannelSheetWidthEnum);
109+ IssmDouble c_t = element->FindParam(HydrologyPressureMeltCoefficientEnum);
110+
111+ Input* h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
112+ Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input);
113+ Input* b_input = element->GetInput(BedEnum); _assert_(b_input);
114+ Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
115+ Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
116+ Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum); _assert_(phiold_input);
117+ Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input);
118+
119+ /*Get tangent vector*/
120+ IssmDouble tx = xyz_list_tria[index2][0] - xyz_list_tria[index1][0];
121+ IssmDouble ty = xyz_list_tria[index2][1] - xyz_list_tria[index1][1];
122+ tx = tx/sqrt(tx*tx+ty*ty);
123+ ty = ty/sqrt(tx*tx+ty*ty);
124+
125+ /*Evaluate fields on center of edge*/
126+ GaussTria* gauss=new GaussTria();
127+ gauss->GaussEdgeCenter(index1,index2);
128+
129+ tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
130+ tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
131+ tria->GetSegmentNodalFunctionsDerivatives(&dbasisdx[0],&xyz_list_tria[0][0],gauss,index1,index2,tria->FiniteElement());
132+ dbasisds[0] = dbasisdx[0*2+0]*tx + dbasisdx[0*2+1]*ty;
133+ dbasisds[1] = dbasisdx[1*2+0]*tx + dbasisdx[1*2+1]*ty;
134+
135+ /*Get input values at gauss points*/
136+ phi_input->GetInputValue(&phi,gauss);
137+ phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
138+ h_input->GetInputValue(&h,gauss);
139+ B_input->GetInputValue(&B,gauss);
140+ n_input->GetInputValue(&n,gauss);
141+ b_input->GetInputValue(&b,gauss);
142+ b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
143+ H_input->GetInputValue(&H,gauss);
144+
145+ /*Get values for a few potentials*/
146+ phi_0 = rho_water*g*b + rho_ice*g*H;
147+ dphids = dphi[0]*tx + dphi[1]*ty;
148+ dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
149+
150+ /*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
151+ IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(fabs(dphids),BETA_C-2.);
152+ IssmDouble Ks = ks * pow(h ,ALPHA_S) * pow(fabs(dphids),BETA_S-2.);
153+
154+ /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
155+ qc = - Ks * dphids;
156+
157+ /*d(phi - phi_m)/ds*/
158+ dPw = dphids - dphimds;
159+
160+ /*Compute f factor*/
161+ fFactor = 0.;
162+ if(this->S>0. || qc*dPw>0.){
163+ fFactor = lc * qc;
164+ }
165+
166+ /*Compute Afactor and Bfactor*/
167+ Afactor = C_W*c_t*rho_water;
168+ Bfactor = 1./L * (1./rho_ice - 1./rho_water);
169+ if(dphids>0){
170+ Xifactor = + Bfactor * (fabs(-Kc*dphids) + fabs(lc*qc));
171+ }
172+ else{
173+ Xifactor = - Bfactor * (fabs(-Kc*dphids) + fabs(lc*qc));
174+ }
175+
176+ _error_("STOP");
177+
178+ /*Diffusive term*/
179+ for(int i=0;i<numnodes;i++){
180+ for(int j=0;j<numnodes;j++){
181+ /*GlaDSCoupledSolver.F90 line 1659*/
182+ Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
183+ +Kc*dbasisds[i]*dbasisds[j] /*Diffusion term*/
184+ - Afactor * Bfactor* Kc * dPw * basis[i] * dbasisds[j] /*First part of Pi*/
185+ + Afactor * fFactor * Bfactor * basis[i] * dbasisds[j] /*Second part of Pi*/
186+ + Xifactor* basis[i] * dbasisds[j] /*Xi term*/
187+ );
188+ }
189+ }
190+
191+ /*Clean up and return*/
192+ delete gauss;
193+}
194+/*}}}*/
195Index: ../trunk-jpl/src/c/classes/Loads/Channel.h
196===================================================================
197--- ../trunk-jpl/src/c/classes/Loads/Channel.h (revision 23964)
198+++ ../trunk-jpl/src/c/classes/Loads/Channel.h (revision 23965)
199@@ -73,7 +73,7 @@
200 void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
201 /*}}}*/
202 /*Channel management:{{{*/
203- void GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
204+ void UpdateChannelCrossSection(void);
205 ElementVector* CreatePVectorHydrologyGlaDS(void);
206 ElementMatrix* CreateKMatrixHydrologyGlaDS(void);
207 /*}}}*/
Note: See TracBrowser for help on using the repository browser.