Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23964) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23965) @@ -372,5 +372,4 @@ void HydrologyGlaDSAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ element->InputUpdateFromSolutionOneDof(solution,HydraulicPotentialEnum); - this->UpdateEffectivePressure(element); }/*}}}*/ void HydrologyGlaDSAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ @@ -510,2 +509,16 @@ xDelete(N); }/*}}}*/ +void HydrologyGlaDSAnalysis::UpdateChannelCrossSection(FemModel* femmodel){/*{{{*/ + + bool ischannels; + femmodel->parameters->FindParam(&ischannels,HydrologyIschannelsEnum); + if(!ischannels) return; + + for(int i=0;iloads->Size();i++){ + if(femmodel->loads->GetEnum(i)==ChannelEnum){ + Channel* channel=(Channel*)femmodel->loads->GetObjectByOffset(i); + channel->UpdateChannelCrossSection(); + } + } + +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h (revision 23964) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h (revision 23965) @@ -34,6 +34,7 @@ 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 Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 23964) +++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 23965) @@ -572,2 +572,125 @@ } /*}}}*/ +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;ivalues[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; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Channel.h (revision 23964) +++ /issm/trunk-jpl/src/c/classes/Loads/Channel.h (revision 23965) @@ -74,5 +74,5 @@ /*}}}*/ /*Channel management:{{{*/ - void GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]); + void UpdateChannelCrossSection(void); ElementVector* CreatePVectorHydrologyGlaDS(void); ElementMatrix* CreateKMatrixHydrologyGlaDS(void); Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 23964) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 23965) @@ -198,6 +198,11 @@ 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; }