Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 17020) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 17021) @@ -247,6 +247,7 @@ scalar = Jdet*gauss->weight*(-transfer); if(dt!=0.) scalar = scalar*dt; - for(int i=0;ivalues[i]+=scalar*basis[i]; - + + for(int i=0;ivalues[i]+=scalar*basis[i]; + /*Transient term*/ if(dt!=0.){ @@ -254,5 +255,6 @@ old_wh_input->GetInputValue(&water_head,gauss); scalar = Jdet*gauss->weight*water_head*epl_specificstoring*epl_thickness; - for(int i=0;ivalues[i]+=scalar*basis[i]; + + for(int i=0;ivalues[i]+=scalar*basis[i]; } } Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 17020) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 17021) @@ -85,4 +85,6 @@ iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); iomodel->FetchDataToInput(elements,SedimentHeadEnum); + iomodel->FetchDataToInput(elements,HydrologydcSedimentTransmitivityEnum); + if(isefficientlayer)iomodel->FetchDataToInput(elements,HydrologydcMaskEplactiveNodeEnum); }/*}}}*/ @@ -168,4 +170,5 @@ /*Intermediaries */ IssmDouble D_scalar,Jdet,dt; + IssmDouble sediment_transmitivity; IssmDouble *xyz_list = NULL; @@ -177,18 +180,18 @@ IssmDouble* B = xNew(2*numnodes); IssmDouble* basis = xNew(numnodes); - IssmDouble D[2][2]={0.}; + IssmDouble D[2][2]= {0.}; /*Retrieve all inputs and parameters*/ - basalelement->GetVerticesCoordinates(&xyz_list); - IssmDouble sediment_storing = SedimentStoring(basalelement); - IssmDouble sediment_transmitivity = basalelement->GetMaterialParameter(HydrologydcSedimentTransmitivityEnum); - basalelement->FindParam(&dt,TimesteppingTimeStepEnum); + basalelement ->GetVerticesCoordinates(&xyz_list); + IssmDouble sediment_storing = SedimentStoring(basalelement); + Input* SedTrans_input=basalelement ->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input); + basalelement ->FindParam(&dt,TimesteppingTimeStepEnum); /* Start looping on the number of gaussian points: */ Gauss* gauss=basalelement->NewGauss(2); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + for(int ig=gauss -> begin();igend();ig++){ + gauss -> GaussPoint(ig); + basalelement -> JacobianDeterminant(&Jdet,xyz_list,gauss); + SedTrans_input -> GetInputValue(&sediment_transmitivity,gauss); /*Diffusivity*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17020) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17021) @@ -4660,11 +4660,10 @@ bool active_element; int transfermethod,step; - IssmDouble sed_trans,sed_thick; IssmDouble leakage,h_max; - IssmDouble wh_trans; + IssmDouble wh_trans,sed_thick; IssmDouble activeEpl[numdof],epl_thickness[numdof]; IssmDouble epl_specificstoring[numdof],sedstoring[numdof]; IssmDouble epl_head[numdof],sed_head[numdof]; - IssmDouble old_transfer[numdof]; + IssmDouble old_transfer[numdof],sed_trans[numdof]; Input* active_element_input=NULL; @@ -4694,4 +4693,5 @@ GetInputListOnVertices(&sed_head[0],SedimentHeadEnum); + GetInputListOnVertices(&sed_trans[0],HydrologydcSedimentTransmitivityEnum); GetInputListOnVertices(&epl_head[0],EplHeadEnum); GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum); @@ -4699,18 +4699,12 @@ this->parameters->FindParam(&leakage,HydrologydcLeakageFactorEnum); - sed_trans = matpar->GetSedimentTransmitivity(); sed_thick = matpar->GetSedimentThickness(); if(!active_element){ /*No transfer if the EPL is not active*/ - /* for(int i=0;iSetValue(doflist[i],wh_trans,INS_VAL); */ - /* } */ } else{ - //GetInputListOnVertices(&old_transfer[0],WaterTransferEnum); + GetInputListOnVertices(&old_transfer[0],WaterTransferEnum); for(int i=0;iGetSedimentStoring(); this->GetHydrologyDCInefficientHmax(&h_max,i); - + + //avoiding transfer at first time + if(epl_head[i]==0.0)epl_head[i]=sed_head[i]; /*EPL head higher than sediment head, transfer from the epl to the sediment*/ if(epl_head[i]>sed_head[i]){ - wh_trans=epl_specificstoring[i]*epl_thickness[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick); - + wh_trans=epl_specificstoring[i]*epl_thickness[i]*sed_trans[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); /*No transfer if the sediment head is allready at the maximum*/ if(sed_head[i]>=h_max){ @@ -4730,12 +4725,14 @@ /*EPL head lower than sediment head, transfer from the sediment to the epl*/ else if(epl_head[i]<=sed_head[i]){ - wh_trans=sedstoring[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick); + wh_trans=sedstoring[i]*sed_trans[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); } /*Introduce relaxation*/ - //wh_trans=old_transfer[i]+0.9*(wh_trans-old_transfer[i]); + wh_trans=old_transfer[i]+0.66*(wh_trans-old_transfer[i]); /*Assign output pointer*/ transfer->SetValue(doflist[i],wh_trans,INS_VAL); + + } } @@ -4881,5 +4878,4 @@ /*And proceed to the real thing*/ thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))); - } } Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 17020) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 17021) @@ -59,5 +59,4 @@ iomodel->Constant(&this->sediment_porosity,HydrologydcSedimentPorosityEnum); iomodel->Constant(&this->sediment_thickness,HydrologydcSedimentThicknessEnum); - iomodel->Constant(&this->sediment_transmitivity,HydrologydcSedimentTransmitivityEnum); iomodel->Constant(&this->water_compressibility,HydrologydcWaterCompressibilityEnum); iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum); @@ -262,5 +261,4 @@ case HydrologydcEplInitialThicknessEnum: return this->epl_init_thickness; case HydrologydcWaterCompressibilityEnum: return this->water_compressibility; - case HydrologydcSedimentTransmitivityEnum: return this->sediment_transmitivity; case HydrologyshreveCREnum: return this->hydro_CR; case HydrologyshreveKnEnum: return this->hydro_kn; @@ -390,10 +388,5 @@ return this->rho_freshwater * this->g * this->epl_porosity * (this->water_compressibility + (this->epl_compressibility / this->epl_porosity)); -} -/*}}}*/ -/*FUNCTION Matpar::GetSedimentTransitivity {{{*/ -IssmDouble Matpar::GetSedimentTransmitivity(){ - return sediment_transmitivity; -} +} /*}}}*/ /*FUNCTION Matpar::GetSedimentThickness {{{*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 17020) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 17021) @@ -44,5 +44,4 @@ IssmDouble sediment_porosity; IssmDouble sediment_thickness; - IssmDouble sediment_transmitivity; IssmDouble water_compressibility; @@ -122,5 +121,4 @@ IssmDouble GetSedimentStoring(); IssmDouble GetEplSpecificStoring(); - IssmDouble GetSedimentTransmitivity(); IssmDouble GetSedimentThickness(); IssmDouble GetEplConductivity(); Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17020) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17021) @@ -97,5 +97,11 @@ sedconverged=false; for(;;){ - + /* if(isefficientlayer){ */ + /* /\*Updating Elemental Mask*\/ */ + /* HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); */ + /* analysis->ElementizeEplMask(femmodel); */ + /* delete analysis; */ + /* femmodel->HydrologyTransferx(); */ + /* } */ uf_sed_sub_iter=uf_sed->Duplicate();_assert_(uf_sed_sub_iter); uf_sed->Copy(uf_sed_sub_iter); @@ -123,11 +129,4 @@ if(sedconverged){ - if(isefficientlayer){ - /*Updating Elemental Mask*/ - HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); - analysis->ElementizeEplMask(femmodel); - delete analysis; - femmodel->HydrologyTransferx(); - } sedconverged=false; /*Checking convegence on the value of the sediment head*/ @@ -149,4 +148,11 @@ if(sedconverged){ + if(isefficientlayer){ + /*Updating Elemental Mask*/ + HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); + analysis->ElementizeEplMask(femmodel); + delete analysis; + femmodel->HydrologyTransferx(); + } femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum); InputUpdateFromConstantx(femmodel,sedconverged,ConvergedEnum); @@ -159,15 +165,19 @@ /*Second layer*/ if(isefficientlayer){ - + femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum); InputUpdateFromConstantx(femmodel,false,ConvergedEnum); femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum); + + /*Iteration on the EPL layer*/ eplconverged = false; for(;;){ - - /*Start by retrieving the EPL head slopes*/ + + //updating mask after the computation of the epl thickness (Allow to close too thin EPL) + femmodel->HydrologyEPLupdateDomainx(); + /*Start by retrieving the EPL head slopes and compute EPL Thickness*/ if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum); @@ -177,15 +187,17 @@ femmodel->parameters->SetParam(EplHeadSlopeYEnum,InputToL2ProjectEnum); solutionsequence_linear(femmodel); - - femmodel->HydrologyEPLThicknessx(); - + femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); - - //updating mask after the computation of the epl thickness (Allow to close too thin EPL) - femmodel->HydrologyEPLupdateDomainx(); + + /*Updating Elemental Mask*/ + HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); + analysis->ElementizeEplMask(femmodel); + delete analysis; + femmodel->HydrologyTransferx(); ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter); ug_epl->Copy(ug_epl_sub_iter); + SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); @@ -213,5 +225,4 @@ eplcount++; - if(eplconverged){ eplconverged=false; @@ -232,10 +243,5 @@ if(eplconverged){ - /*Updating Elemental Mask*/ - HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); - analysis->ElementizeEplMask(femmodel); - delete analysis; - femmodel->HydrologyTransferx(); - + InputUpdateFromConstantx(femmodel,eplconverged,ConvergedEnum); InputUpdateFromConstantx(femmodel,sediment_kmax,MeltingOffsetEnum);