Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 22840) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 22841) @@ -74,5 +74,5 @@ bool isefficientlayer; int hydrology_model; - + /*Fetch data needed: */ iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); @@ -161,5 +161,5 @@ loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum)); loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum)); - } + } } } @@ -233,5 +233,5 @@ active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); } - + /* Start looping on the number of gaussian points: */ Gauss* gauss=basalelement->NewGauss(2); @@ -250,5 +250,5 @@ D[0][0]=D_scalar; D[1][1]=D_scalar; - GetB(B,basalelement,xyz_list,gauss); + GetB(B,basalelement,xyz_list,gauss); TripleMultiply(B,2,numnodes,1, &D[0][0],2,2,0, @@ -265,5 +265,5 @@ basis,1,numnodes,0, &Ke->values[0],1); - + /*Transfer EPL part*/ if(isefficientlayer){ @@ -288,5 +288,5 @@ delete gauss; if(domaintype!=Domain2DhorizontalEnum){ - basalelement->DeleteMaterials(); + basalelement->DeleteMaterials(); delete basalelement; } @@ -412,5 +412,5 @@ delete gauss; if(domaintype!=Domain2DhorizontalEnum){ - basalelement->DeleteMaterials(); + basalelement->DeleteMaterials(); delete basalelement; } @@ -419,7 +419,7 @@ void HydrologyDCInefficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. + /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. * For node i, Bi can be expressed in the actual coordinate system - * by: + * by: * Bi=[ dN/dx ] * [ dN/dy ] @@ -507,5 +507,5 @@ IssmDouble rho_ice = basalelement->GetMaterialParameter(MaterialsRhoIceEnum); IssmDouble g = basalelement->GetMaterialParameter(ConstantsGEnum); - + basalelement->GetInputListOnVertices(thickness,ThicknessEnum); basalelement->GetInputListOnVertices(base,BaseEnum); @@ -542,5 +542,5 @@ xDelete(doflist); if(domaintype!=Domain2DhorizontalEnum){ - basalelement->DeleteMaterials(); + basalelement->DeleteMaterials(); delete basalelement; } @@ -554,5 +554,5 @@ IssmDouble HydrologyDCInefficientAnalysis::SedimentStoring(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input){/*{{{*/ int unconf_scheme; - IssmDouble expfac; + IssmDouble expfac; IssmDouble sediment_storing; IssmDouble storing,yield; @@ -610,5 +610,5 @@ sed_head_input->GetInputValue(&prestep_head,gauss); water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness))); - + //min definition of the if test sediment_transmitivity=FullLayer_transmitivity*min(water_sheet/sediment_thickness,1.); @@ -625,5 +625,5 @@ void HydrologyDCInefficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/ - + int hmax_flag; IssmDouble h_max; @@ -632,5 +632,5 @@ /*Get the flag to the limitation method*/ element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); - + /*Switch between the different cases*/ switch(hmax_flag){ @@ -674,5 +674,5 @@ case 1: element->FindParam(&leakage,HydrologydcLeakageFactorEnum); - transfer=+leakage; + transfer=+leakage; break; default: @@ -698,7 +698,4 @@ _assert_(epl_head_input); epl_head_input->GetInputValue(&epl_head,gauss); - if (element->Id()==42){ - _printf_("epl head in sed Pvec transfer is "<< epl_head <<"\n"); - } element->FindParam(&leakage,HydrologydcLeakageFactorEnum); transfer=+epl_head*leakage; @@ -718,5 +715,5 @@ element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input); - + if(node_mask_input->Max()>0.){ element_active = true; Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 22840) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 22841) @@ -13,5 +13,5 @@ #include "../classes.h" #include "shared/shared.h" -/*}}}*/ +/*}}}*/ /*Constructors/destructors*/ @@ -88,4 +88,5 @@ switch(CoupledFlag){ + //This should be removed at some point and the Enum uniformalized case 0: Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel)); @@ -95,5 +96,5 @@ break; case 2: - element->GetInputValue(&Neff,gauss,EffectivePressureEnum); + element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum); break; default: @@ -174,6 +175,6 @@ void Friction::GetAlphaViscousComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/ - /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p. - * FrictionGetAlphaComplement is used in control methods on drag, and it computes: + /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p. + * FrictionGetAlphaComplement is used in control methods on drag, and it computes: * alpha_complement= Neff ^r * vel ^s*/ @@ -207,10 +208,10 @@ case 0: Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel)); - break; + break; case 1: element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum); break; case 2: - element->GetInputValue(&Neff,gauss,EffectivePressureEnum); + element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum); break; default: @@ -285,5 +286,5 @@ void Friction::GetAlpha2Coulomb(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ - /*This routine calculates the basal friction coefficient + /*This routine calculates the basal friction coefficient alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/ @@ -319,10 +320,10 @@ case 0: Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel)); - break; + break; case 1: element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum); break; case 2: - element->GetInputValue(&Neff,gauss,EffectivePressureEnum); + element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum); break; default: @@ -370,5 +371,5 @@ void Friction::GetAlpha2Hydro(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ - /*This routine calculates the basal friction coefficient + /*This routine calculates the basal friction coefficient Based on Gagliardini 2007, needs a good effective pressure computation Not tested so far so use at your own risks @@ -405,14 +406,14 @@ IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum); - + switch(CoupledFlag){ case 0: Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel)); - break; + break; case 1: element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum); break; case 2: - element->GetInputValue(&Neff,gauss,EffectivePressureEnum); + element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum); break; default: @@ -453,5 +454,5 @@ Gamma=(Chi/(1. + alpha * pow(Chi,q_exp))); /*Check to prevent dividing by zero if vmag==0*/ - if(vmag==0.) alpha2=0.; + if(vmag==0.) alpha2=0.; else if (Neff==0) alpha2=0.0; else alpha2=Neff * C_param * pow(Gamma,1./n) * 1/vmag; @@ -570,5 +571,5 @@ void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ - /*This routine calculates the basal friction coefficient + /*This routine calculates the basal friction coefficient alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/ @@ -607,5 +608,5 @@ break; case 2: - element->GetInputValue(&Neff,gauss,EffectivePressureEnum); + element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum); break; default: @@ -644,5 +645,5 @@ void Friction::GetAlpha2WaterLayer(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ - /*This routine calculates the basal friction coefficient + /*This routine calculates the basal friction coefficient alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/ Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 22840) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 22841) @@ -1,5 +1,5 @@ /*!\file: hydrology_core.cpp - * \brief: core of the hydrology solution - */ + * \brief: core of the hydrology solution + */ #include "./cores.h" @@ -23,5 +23,5 @@ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum); - femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); femmodel->parameters->FindParam(&numoutputs,HydrologyNumRequestedOutputsEnum); if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,HydrologyRequestedOutputsEnum); @@ -36,8 +36,8 @@ femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum); solutionsequence_nonlinear(femmodel,modify_loads); - + /*transfer water column thickness to old water column thickness: */ InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum); - + } @@ -90,6 +90,4 @@ /*save preceding timestep*/ InputDuplicatex(femmodel,SedimentHeadEnum,SedimentHeadOldEnum); - InputDuplicatex(femmodel,EplHeadEnum,EplHeadOldEnum); - InputDuplicatex(femmodel,HydrologydcEplThicknessEnum,HydrologydcEplThicknessOldEnum); /*Proceed now to heads computations*/ solutionsequence_hydro_nonlinear(femmodel); @@ -110,8 +108,8 @@ } } - else if (hydrology_model==HydrologysommersEnum){ - femmodel->SetCurrentConfiguration(HydrologySommersAnalysisEnum); - InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum); - solutionsequence_shakti_nonlinear(femmodel); + else if (hydrology_model==HydrologysommersEnum){ + femmodel->SetCurrentConfiguration(HydrologySommersAnalysisEnum); + InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum); + solutionsequence_shakti_nonlinear(femmodel); if(VerboseSolution()) _printf0_(" updating gap height\n"); HydrologySommersAnalysis* analysis = new HydrologySommersAnalysis(); @@ -126,5 +124,5 @@ femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); } - /*Free ressources:*/ + /*Free ressources:*/ if(numoutputs){ for(int i=0;i