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

../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp

@@ -371 +371 @@
 }/*}}}*/
 void           HydrologyGlaDSAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
         element->InputUpdateFromSolutionOneDof(solution,HydraulicPotentialEnum);
-        this->UpdateEffectivePressure(element);
 }/*}}}*/
 void           HydrologyGlaDSAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
         /*Default, do nothing*/
@@ -509 +508 @@
         delete gauss;
         xDelete<IssmDouble>(N);
 }/*}}}*/
+void HydrologyGlaDSAnalysis::UpdateChannelCrossSection(FemModel* femmodel){/*{{{*/
+
+        bool ischannels;
+        femmodel->parameters->FindParam(&ischannels,HydrologyIschannelsEnum);
+        if(!ischannels) return;
+
+        for(int i=0;i<femmodel->loads->Size();i++){
+                if(femmodel->loads->GetEnum(i)==ChannelEnum){
+                        Channel* channel=(Channel*)femmodel->loads->GetObjectByOffset(i);
+                        channel->UpdateChannelCrossSection();
+                }
+        }
+
+}/*}}}*/

../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h

@@ -33 +33 @@
                 /*Specific to GlaDS*/
                 void UpdateSheetThickness(FemModel* femmodel);
                 void UpdateSheetThickness(Element*  element);
+                void UpdateChannelCrossSection(FemModel* femmodel);
                 void UpdateEffectivePressure(FemModel* femmodel);
-                void UpdateEffectivePressure(Element*  element);
+                void UpdateEffectivePressure(Element* element);
 };
 #endif

../trunk-jpl/src/c/cores/hydrology_core.cpp

@@ -197 +197 @@
                 InputDuplicatex(femmodel,HydraulicPotentialEnum,HydraulicPotentialOldEnum);
                 analysis->UpdateEffectivePressure(femmodel);
                 solutionsequence_shakti_nonlinear(femmodel);
-                if(VerboseSolution()) _printf0_("   updating sheet thickness\n");
+
+                if(VerboseSolution()) _printf0_("   updating effective pressure\n");
+                analysis->UpdateEffectivePressure(femmodel);
+                if(VerboseSolution()) _printf0_("   updating sheet thickness\n"); /*Uses N, so needs to come after*/
                 analysis->UpdateSheetThickness(femmodel);
+                if(VerboseSolution()) _printf0_("   updating channels cross section\n");
+                analysis->UpdateChannelCrossSection(femmodel);
                 delete analysis;
         }
 

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

@@ -571 +571 @@
         return pe;
 }
 /*}}}*/
+void           Channel::UpdateChannelCrossSection(void){/*{{{*/
+
+        /*Initialize Element matrix and return if necessary*/
+        Tria*  tria=(Tria*)element;
+        if(!tria->IsIceInElement()){
+                this->S = 0.;
+                return;
+        }
+        _assert_(tria->FiniteElement()==P1Enum); 
+        int index1=tria->GetNodeIndex(nodes[0]);
+        int index2=tria->GetNodeIndex(nodes[1]);
+
+        /*Intermediaries */
+        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;
+        IssmDouble  xyz_list[NUMVERTICES][3];
+        IssmDouble  xyz_list_tria[3][3];
+        const int   numnodes = NUMNODES;
+
+        /*Initialize Element vector and other vectors*/
+        ElementMatrix* Ke=new ElementMatrix(this->nodes,NUMNODES,this->parameters);
+        IssmDouble     basis[NUMNODES];
+        IssmDouble     dbasisdx[2*NUMNODES];
+        IssmDouble     dbasisds[NUMNODES];
+
+        /*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 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);
+        Input* b_input      = element->GetInput(BedEnum); _assert_(b_input);
+        Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
+        Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
+        Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
+        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);
+
+        /*Evaluate fields on center of edge*/
+        GaussTria* gauss=new GaussTria();
+        gauss->GaussEdgeCenter(index1,index2);
+
+        tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
+        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);
+        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;
+        }
+
+        /*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));
+        }
+
+        _error_("STOP");
+
+        /*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]                               /*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*/
+                                                );
+                }
+        }
+
+        /*Clean up and return*/
+        delete gauss;
+}
+/*}}}*/

../trunk-jpl/src/c/classes/Loads/Channel.h

@@ -73 +73 @@
                 void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
                 /*}}}*/
                 /*Channel management:{{{*/
-                void           GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
+                void           UpdateChannelCrossSection(void);
                 ElementVector* CreatePVectorHydrologyGlaDS(void);
                 ElementMatrix* CreateKMatrixHydrologyGlaDS(void);
                 /*}}}*/