Changeset 23963

Timestamp:

05/30/19 14:53:00 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on channels

File:

• issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/classes/Loads/Channel.cpp

@@ -14 +14 @@
 #include "../classes.h"
 /*}}}*/ 
+
+/*Macros*/
 #define NUMNODES    2
 #define NUMVERTICES 2
+
+#define C_W         4.22e3   /*specific heat capacity of water (J/kg/K)*/
+#define ALPHA_C     5./4.
+#define BETA_C      3./2.
+/*Make sure these are the same as in HydrologyGlaDSAnalysis::CreateKMatrix*/
+#define ALPHA_S     5./4.
+#define BETA_S      3./2.
 
 /*Channel constructors and destructor*/
@@ -338 +347 @@
 
         /*Intermediaries */
-        IssmDouble  Jdet,v1;
-        IssmDouble  A,B,n,phi_old,phi,phi_0;
-        IssmDouble  H,h,b;
+        IssmDouble  Jdet,v1,qc,fFactor;
+        IssmDouble  A,B,n,phi_old,phi,phi_0,dPw;
+        IssmDouble  H,h,b,dphi[2],dphids,dphimds,db[2],dbds;
         IssmDouble  xyz_list[NUMVERTICES][3];
         IssmDouble  xyz_list_tria[3][3];
@@ -354 +363 @@
         GetVerticesCoordinates(&xyz_list[0][0]     ,this->vertices,NUMVERTICES);
         GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3);
+
         IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
@@ -359 +369 @@
         IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
         IssmDouble g         = element->FindParam(ConstantsGEnum);
+        IssmDouble kc        = element->FindParam(HydrologyChannelConductivityEnum);
+        IssmDouble ks        = element->FindParam(HydrologySheetConductivityEnum);
+        IssmDouble lc        = element->FindParam(HydrologyChannelSheetWidthEnum);
+        IssmDouble c_t       = element->FindParam(HydrologyPressureMeltCoefficientEnum);
+
         Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
         Input* H_input      = element->GetInput(ThicknessEnum); _assert_(H_input);
@@ -381 +396 @@
                 tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
                 tria->GetSegmentNodalFunctionsDerivatives(&dbasisdx[0],&xyz_list_tria[0][0],gauss,index1,index2,tria->FiniteElement());
+                dbasisds[0] = dbasisdx[0*2+0]*tx + dbasisdx[0*2+1]*ty;
+                dbasisds[1] = dbasisdx[1*2+0]*tx + dbasisdx[1*2+1]*ty;
 
                 /*Get input values at gauss points*/
+                phi_input->GetInputValue(&phi,gauss);
+                phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
                 h_input->GetInputValue(&h,gauss);
                 B_input->GetInputValue(&B,gauss);
                 n_input->GetInputValue(&n,gauss);
-                phi_input->GetInputValue(&phi,gauss);
                 b_input->GetInputValue(&b,gauss);
+                b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
                 H_input->GetInputValue(&H,gauss);
 
+                /*Get values for a few potentials*/
+                phi_0   = rho_water*g*b + rho_ice*g*H;
+                dphids  = dphi[0]*tx + dphi[1]*ty;
+                dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
+
+                /*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
+                IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(fabs(dphids),BETA_C-2.);
+                IssmDouble Ks = ks * pow(h      ,ALPHA_S) * pow(fabs(dphids),BETA_S-2.);
+
+                /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
+                qc = - Ks * dphids;
+
+                /*d(phi - phi_m)/ds*/
+                dPw = dphids - dphimds;
+
+                /*Compute f factor*/
+                fFactor = 0.;
+                if(this->S>0. || qc*dPw>0.){
+                        fFactor = lc * qc;
+                }
+
+                /*Diffusive term*/
+                for(int i=0;i<numnodes;i++){
+                        for(int j=0;j<numnodes;j++){
+                                Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
+                                                        +Kc*dbasisds[i]*dbasisds[j] /*GlaDSCoupledSolver.F90 line 1659*/
+                                                        );
+                        }
+                }
+
                 /*Closing rate term, see Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/
-                phi_0 = rho_water*g*b + rho_ice*g*H;
                 A=pow(B,-n);
                 v1 = 2./pow(n,n)*A*S*(pow(fabs(phi_0 - phi),n-1.)*( - n));