Changeset 23964

Timestamp:

05/30/19 15:30:59 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: finished RHS

File:

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

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

@@ -347 +347 @@
 
         /*Intermediaries */
-        IssmDouble  Jdet,v1,qc,fFactor;
+        IssmDouble  Jdet,v1,qc,fFactor,Afactor,Bfactor,Xifactor;
         IssmDouble  A,B,n,phi_old,phi,phi_0,dPw;
         IssmDouble  H,h,b,dphi[2],dphids,dphimds,db[2],dbds;
@@ -367 +367 @@
         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
         IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
-        IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
         IssmDouble g         = element->FindParam(ConstantsGEnum);
         IssmDouble kc        = element->FindParam(HydrologyChannelConductivityEnum);
@@ -430 +429 @@
                 }
 
+                /*Compute Afactor and Bfactor*/
+                Afactor = C_W*c_t*rho_water;
+                Bfactor = 1./L * (1./rho_ice - 1./rho_water);
+                if(dphids>0){
+                        Xifactor = + Bfactor * (fabs(-Kc*dphids) + fabs(lc*qc));
+                }
+                else{
+                        Xifactor = - Bfactor * (fabs(-Kc*dphids) + fabs(lc*qc));
+                }
+
                 /*Diffusive term*/
                 for(int i=0;i<numnodes;i++){
                         for(int j=0;j<numnodes;j++){
+                                /*GlaDSCoupledSolver.F90 line 1659*/
                                 Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
-                                                        +Kc*dbasisds[i]*dbasisds[j] /*GlaDSCoupledSolver.F90 line 1659*/
+                                                        +Kc*dbasisds[i]*dbasisds[j]                               /*Diffusion term*/
+                                                        - Afactor * Bfactor* Kc * dPw * basis[i] * dbasisds[j]    /*First part of Pi*/
+                                                        + Afactor * fFactor * Bfactor * basis[i] * dbasisds[j]    /*Second part of Pi*/
+                                                        + Xifactor* basis[i] * dbasisds[j]                        /*Xi term*/
                                                         );
                         }
@@ -447 +460 @@
                         }
                 }
-
-                /*Transient term if dt>0*/
-                if(dt>0.){
-                }
         }
 
@@ -468 +477 @@
 
         /*Intermediaries */
-        IssmDouble  Jdet,v2;
-        IssmDouble  A,B,n,phi_old,phi,phi_0;
-        IssmDouble  H,h,b;
+        IssmDouble  Jdet,v2,Afactor,Bfactor,fFactor;
+        IssmDouble  A,B,n,phi_old,phi,phi_0,dphimds;
+        IssmDouble  H,h,b,db[2],dphids,qc,dPw;
         IssmDouble  xyz_list[NUMVERTICES][3];
+        IssmDouble  xyz_list_tria[3][3];
         const int   numnodes = NUMNODES;
 
@@ -480 +490 @@
         /*Retrieve all inputs and parameters*/
         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);
         IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
-        IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
+        IssmDouble kc        = element->FindParam(HydrologyChannelConductivityEnum);
+        IssmDouble ks        = element->FindParam(HydrologySheetConductivityEnum);
         IssmDouble g         = element->FindParam(ConstantsGEnum);
+        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);
@@ -493 +509 @@
         Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
 
+        /*Get tangent vector*/
+        IssmDouble tx = xyz_list_tria[index2][0] - xyz_list_tria[index1][0];
+        IssmDouble ty = xyz_list_tria[index2][1] - xyz_list_tria[index1][1];
+        tx = tx/sqrt(tx*tx+ty*ty);
+        ty = ty/sqrt(tx*tx+ty*ty);
+
         /* Start  looping on the number of gaussian points: */
         Gauss* gauss=new GaussTria(index1,index2,2);
@@ -502 +524 @@
 
                 /*Get input values at gauss points*/
+                b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
                 h_input->GetInputValue(&h,gauss);
                 B_input->GetInputValue(&B,gauss);
@@ -509 +532 @@
                 H_input->GetInputValue(&H,gauss);
 
+                /*Get values for a few potentials*/
+                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;
+                }
+
+                /*Compute Afactor and Bfactor*/
+                Afactor = C_W*c_t*rho_water;
+                Bfactor = 1./L * (1./rho_ice - 1./rho_water);
+
                 /*Compute closing rate*/
                 /*See Gagliardini and Werder 2018 eq. A2 (v = v2(phi_i) + v1*phi_{i+1})*/
@@ -515 +561 @@
                 v2 = 2./pow(n,n)*A*this->S*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi));
 
-                for(int i=0;i<numnodes;i++) pe->values[i]+= - Jdet*gauss->weight*(-v2)*basis[i];
-
-                /*Transient term if dt>0*/
-                if(dt>0.){
-                        //phiold_input->GetInputValue(&phi_old,gauss);
-                        //for(int i=0;i<numnodes;i++) pe->values[i] += gauss->weight*Jdet*e_v/(rho_water*g*dt)*phi_old*basis[i];
+                for(int i=0;i<numnodes;i++){
+                        pe->values[i]+= - Jdet*gauss->weight*(-v2)*basis[i];
+                        pe->values[i]+= + Jdet*gauss->weight*Afactor*Bfactor*fFactor*dphimds*basis[i];
                 }
         }