Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23970) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23971) @@ -122,5 +122,6 @@ iomodel->FetchDataToInput(elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum); iomodel->FetchDataToInput(elements,"md.hydrology.sheet_conductivity",HydrologySheetConductivityEnum); - iomodel->FetchDataToInput(elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum,0.); + iomodel->FetchDataToInput(elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum); + iomodel->FetchDataToInput(elements,"md.hydrology.moulin_input",HydrologyMoulinInputEnum); iomodel->FetchDataToInput(elements,"md.initialization.watercolumn",HydrologySheetThicknessEnum); iomodel->FetchDataToInput(elements,"md.initialization.hydraulic_potential",HydraulicPotentialEnum); Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 23970) +++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 23971) @@ -348,5 +348,5 @@ /*Intermediaries */ IssmDouble Jdet,v1,qc,fFactor,Afactor,Bfactor,Xifactor; - IssmDouble A,B,n,phi_old,phi,phi_0,dPw; + IssmDouble A,B,n,phi_old,phi,phi_0,dPw,ks,Ngrad; IssmDouble H,h,b,dphi[2],dphids,dphimds,db[2],dbds; IssmDouble xyz_list[NUMVERTICES][3]; @@ -369,5 +369,4 @@ 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); @@ -378,4 +377,5 @@ Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* ks_input = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input); Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum); _assert_(phiold_input); Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); @@ -402,4 +402,5 @@ phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss); h_input->GetInputValue(&h,gauss); + ks_input->GetInputValue(&ks,gauss); B_input->GetInputValue(&B,gauss); n_input->GetInputValue(&n,gauss); @@ -412,8 +413,10 @@ dphids = dphi[0]*tx + dphi[1]*ty; dphimds = rho_water*g*(db[0]*tx + db[1]*ty); + Ngrad = fabs(dphids); + if(Ngrad<1.e-12) Ngrad = 1.e-12; /*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.); + IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(Ngrad,BETA_C-2.); + IssmDouble Ks = ks * pow(h ,ALPHA_S) * pow(Ngrad,BETA_S-2.); /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/ @@ -479,5 +482,5 @@ 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 H,h,b,db[2],dphids,qc,dPw,ks; IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble xyz_list_tria[3][3]; @@ -496,5 +499,4 @@ IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); IssmDouble kc = element->FindParam(HydrologyChannelConductivityEnum); - IssmDouble ks = element->FindParam(HydrologySheetConductivityEnum); IssmDouble g = element->FindParam(ConstantsGEnum); IssmDouble lc = element->FindParam(HydrologyChannelSheetWidthEnum); @@ -506,4 +508,5 @@ Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* ks_input = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input); Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum); _assert_(phiold_input); Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); @@ -526,4 +529,5 @@ b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss); h_input->GetInputValue(&h,gauss); + ks_input->GetInputValue(&ks,gauss); B_input->GetInputValue(&B,gauss); n_input->GetInputValue(&n,gauss); @@ -585,5 +589,5 @@ /*Intermediaries */ - IssmDouble A,B,n,phi,phi_0; + IssmDouble A,B,n,phi,phi_0,ks,Ngrad; IssmDouble H,h,b,dphi[2],dphids,dphimds,db[2],dbds; IssmDouble xyz_list[NUMVERTICES][3]; @@ -599,5 +603,4 @@ 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); @@ -609,4 +612,5 @@ Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* ks_input = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input); Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum); _assert_(phiold_input); Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); @@ -626,4 +630,5 @@ phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss); h_input->GetInputValue(&h,gauss); + ks_input->GetInputValue(&ks,gauss); B_input->GetInputValue(&B,gauss); n_input->GetInputValue(&n,gauss); @@ -636,8 +641,10 @@ dphids = dphi[0]*tx + dphi[1]*ty; dphimds = rho_water*g*(db[0]*tx + db[1]*ty); + Ngrad = fabs(dphids); + if(Ngrad<1.e-12) Ngrad = 1.e-12; /*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.); + IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(Ngrad,BETA_C-2.); + IssmDouble Ks = ks * pow(h ,ALPHA_S) * pow(Ngrad,BETA_S-2.); /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23970) +++ /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23971) @@ -159,4 +159,7 @@ switch(analysis_type){ + case HydrologyGlaDSAnalysisEnum: + pe = this->CreatePVectorHydrologyGlaDS(); + break; case HydrologyShaktiAnalysisEnum: pe = this->CreatePVectorHydrologyShakti(); @@ -310,4 +313,22 @@ /*}}}*/ +ElementVector* Moulin::CreatePVectorHydrologyGlaDS(void){/*{{{*/ + + /*If this node is not the master node (belongs to another partition of the + * mesh), don't add the moulin input a second time*/ + if(node->IsClone()) return NULL; + + IssmDouble moulin_load; + + /*Initialize Element matrix*/ + ElementVector* pe=new ElementVector(&node,1,this->parameters); + + this->element->GetInputValue(&moulin_load,node,HydrologyMoulinInputEnum); + pe->values[0]=moulin_load; + + /*Clean up and return*/ + return pe; +} +/*}}}*/ ElementVector* Moulin::CreatePVectorHydrologyShakti(void){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Moulin.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Moulin.h (revision 23970) +++ /issm/trunk-jpl/src/c/classes/Loads/Moulin.h (revision 23971) @@ -77,4 +77,5 @@ ElementVector* CreatePVectorHydrologyShakti(void); + ElementVector* CreatePVectorHydrologyGlaDS(void); ElementVector* CreatePVectorHydrologyDCInefficient(void); ElementVector* CreatePVectorHydrologyDCEfficient(void);