Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 16953) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 16954) @@ -392,7 +392,10 @@ int numvertices = element->GetNumberOfVertices(); Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); - if(node_mask_input->Max()>0.) element_active = true; - else element_active = false; - + if(node_mask_input->Max()>0.) { + element_active = true; + } + else{ + element_active = false; + } element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active)); } Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16953) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16954) @@ -276,6 +276,6 @@ #ifdef _HAVE_HYDROLOGY_ + virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index)=0; virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0; - virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, int index)=0; virtual void GetHydrologyTransfer(Vector* transfer)=0; virtual void HydrologyEPLGetMask(Vector* mask)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16953) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16954) @@ -6544,5 +6544,4 @@ transfer_input->GetInputValue(&transfer,gauss); scalar = Jdet*gauss->weight*(water_load+transfer); - //printf("are we loading: load,%e ,transfer,%e\n",water_load, transfer); if(reCast(dt)) scalar = scalar*dt; for(int i=0;ivalues[i]+=scalar*basis[i]; @@ -6726,41 +6725,6 @@ /*}}}*/ /*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/ -void Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){ - - int hmax_flag; - IssmDouble h_max; - IssmDouble rho_ice,rho_water; - IssmDouble thickness,bed; - /*Get the flag to the limitation method*/ - this->parameters->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); - - /*Switch between the different cases*/ - switch(hmax_flag){ - case 0: - h_max=1.0e+10; - break; - case 1: - parameters->FindParam(&h_max,HydrologydcSedimentlimitEnum); - break; - case 2: - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoFreshwater(); - this->GetInputValue(&thickness,innode,ThicknessEnum); - this->GetInputValue(&bed,innode,BedEnum); - h_max=((rho_ice*thickness)/rho_water)+bed; - break; - case 3: - _error_("Using normal stress not supported yet"); - break; - default: - _error_("no case higher than 3 for SedimentlimitFlag"); - } - /*Assign output pointer*/ - *ph_max=h_max; -} -/*}}}*/ -/*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/ void Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){ - + int hmax_flag; IssmDouble h_max; @@ -6769,25 +6733,60 @@ /*Get the flag to the limitation method*/ this->parameters->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); - + /*Switch between the different cases*/ switch(hmax_flag){ - case 0: - h_max=1.0e+10; - break; - case 1: - parameters->FindParam(&h_max,HydrologydcSedimentlimitEnum); - break; - case 2: - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoFreshwater(); - this->GetInputValue(&thickness,this->nodes[index],ThicknessEnum); - this->GetInputValue(&bed,this->nodes[index],BedEnum); - h_max=((rho_ice*thickness)/rho_water)+bed; - break; - case 3: - _error_("Using normal stress not supported yet"); - break; - default: - _error_("no case higher than 3 for SedimentlimitFlag"); + case 0: + h_max=1.0e+10; + break; + case 1: + parameters->FindParam(&h_max,HydrologydcSedimentlimitEnum); + break; + case 2: + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoFreshwater(); + this->GetInputValue(&thickness,this->nodes[index],ThicknessEnum); + this->GetInputValue(&bed,this->nodes[index],BedEnum); + h_max=((rho_ice*thickness)/rho_water)+bed; + break; + case 3: + _error_("Using normal stress not supported yet"); + break; + default: + _error_("no case higher than 3 for SedimentlimitFlag"); + } + /*Assign output pointer*/ + *ph_max=h_max; +} +/*}}}*/ +/*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/ +void Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){ + + int hmax_flag; + IssmDouble h_max; + IssmDouble rho_ice,rho_water; + IssmDouble thickness,bed; + /*Get the flag to the limitation method*/ + this->parameters->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); + + /*Switch between the different cases*/ + switch(hmax_flag){ + case 0: + h_max=1.0e+10; + break; + case 1: + parameters->FindParam(&h_max,HydrologydcSedimentlimitEnum); + break; + case 2: + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoFreshwater(); + this->GetInputValue(&thickness,innode,ThicknessEnum); + this->GetInputValue(&bed,innode,BedEnum); + h_max=((rho_ice*thickness)/rho_water)+bed; + break; + case 3: + _error_("Using normal stress not supported yet"); + break; + default: + _error_("no case higher than 3 for SedimentlimitFlag"); } /*Assign output pointer*/ @@ -6831,5 +6830,4 @@ active_element_input->GetInputValue(&active_element); - GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveNodeEnum); GetInputListOnVertices(&sed_head[0],SedimentHeadEnum); GetInputListOnVertices(&epl_head[0],EplHeadEnum); @@ -6861,5 +6859,5 @@ /*No transfer if the sediment head is allready at the maximum*/ - this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]); + this->GetHydrologyDCInefficientHmax(&h_max,i); if(sed_head[i]>=h_max)wh_trans=0.0; } @@ -6870,4 +6868,9 @@ /*Assign output pointer*/ transfer->SetValue(doflist[i],wh_trans,INS_VAL); + /* if(nodes[i]->id>=54){ */ + /* printf("%i %e %e %e \n",nodes[i]->id-54,wh_trans,sed_head[i],epl_head[i]); */ + /* } */ + /* else{*/ + /* printf("%i %e %e %e \n",nodes[i]->id,wh_trans,sed_head[i],epl_head[i]); */ } } @@ -6890,5 +6893,5 @@ GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum); - + for(int i=0;iSetValue(nodes[i]->Sid(),0.,INS_VAL); - } + /* else if(epl_thickness[i]<0.0){ */ + /* vec_mask->SetValue(nodes[i]->Sid(),0.,INS_VAL); */ + /* } */ /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/ - this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]); + this->GetHydrologyDCInefficientHmax(&h_max,i); if(eplhead[i]>=h_max && this->AnyActive()){ for(j=0;j* transfer); void HydrologyEPLGetActive(Vector* active_vec); Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 16953) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 16954) @@ -94,6 +94,6 @@ if(VerboseSolution()) _printf0_(" saving results \n"); if(isefficientlayer){ - int outputs[8] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveNodeEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum,HydrologydcEplThicknessEnum,WaterTransferEnum}; - femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],8); + int outputs[9] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveNodeEnum,HydrologydcMaskEplactiveEltEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum,HydrologydcEplThicknessEnum,WaterTransferEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],9); } else{ Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 16953) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 16954) @@ -20,8 +20,13 @@ Vector* uf_epl=NULL; Vector* uf_epl_sub_iter=NULL; + Vector* ug_epl_sub_iter=NULL; Vector* ug_epl_main_iter=NULL; + Vector* ys=NULL; Vector* dug=NULL; + + //testing stuff + Vector* duf=NULL; Matrix* Kff=NULL; @@ -52,9 +57,16 @@ /*Retrieve inputs as the initial state for the non linear iteration*/ - GetSolutionFromInputsx(&ug_sed,femmodel); + GetSolutionFromInputsx(&ug_sed,femmodel); + + //test + GetSolutionFromInputsx(&uf_sed,femmodel);_assert_(uf_sed); if(isefficientlayer) { femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); GetSolutionFromInputsx(&ug_epl,femmodel); + + //test + GetSolutionFromInputsx(&uf_epl,femmodel);_assert_(uf_epl); + /*Initialize the transfer input*/ HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); @@ -72,7 +84,16 @@ ug_sed_main_iter=ug_sed->Duplicate(); ug_sed->Copy(ug_sed_main_iter); + + //test + uf_sed_sub_iter=uf_sed->Duplicate(); + uf_sed->Copy(uf_sed_sub_iter); + if(isefficientlayer){ ug_epl_main_iter=ug_epl->Duplicate(); ug_epl->Copy(ug_epl_main_iter); + //test + ug_epl_sub_iter=ug_epl->Duplicate(); + ug_epl->Copy(ug_epl_sub_iter); + } @@ -95,7 +116,7 @@ Solverx(&uf_sed,Kff,pf,uf_sed_sub_iter,df,femmodel->parameters); delete Kff; delete pf; delete df; - delete uf_sed_sub_iter; - uf_sed_sub_iter=uf_sed->Duplicate(); - uf_sed->Copy(uf_sed_sub_iter); + /* delete uf_sed_sub_iter; */ + /* uf_sed_sub_iter=uf_sed->Duplicate(); */ + /* uf_sed->Copy(uf_sed_sub_iter); */ delete ug_sed; Mergesolutionfromftogx(&ug_sed,uf_sed,ys,femmodel->nodes,femmodel->parameters); delete ys; @@ -110,5 +131,26 @@ } } + sedcount++; + + //testing stuff + if(sedconverged){ + sedconverged=false; + duf=uf_sed_sub_iter->Duplicate();_assert_(duf); + uf_sed_sub_iter->Copy(duf);_assert_(uf_sed_sub_iter); + duf->AYPX(uf_sed,-1.0); + ndu_sed=duf->Norm(NORM_TWO); + delete duf; + nu_sed=uf_sed_sub_iter->Norm(NORM_TWO); + if (xIsNan(ndu_sed) || xIsNan(nu_sed)) _error_("convergence criterion is NaN!"); + if((ndu_sed/nu_sed)Duplicate();_assert_(uf_sed_sub_iter); + uf_sed->Copy(uf_sed_sub_iter);_assert_(uf_sed); + //end of the crap if(sedconverged){ @@ -163,5 +205,5 @@ Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters); delete Kff; delete pf; delete df; - delete uf_epl_sub_iter; + delete uf_epl_sub_iter; uf_epl_sub_iter=uf_epl->Duplicate(); uf_epl->Copy(uf_epl_sub_iter); @@ -171,4 +213,9 @@ ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); femmodel->HydrologyEPLupdateDomainx(); + /* /\*Updating Nodal Mask*\/ */ + /* HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); */ + /* analysis->ElementizeEplMask(femmodel); */ + /* delete analysis; */ + /* femmodel->HydrologyTransferx(); */ if (!eplconverged){ @@ -182,4 +229,23 @@ } eplcount++; + + //testing stuff + if(eplconverged){ + eplconverged=false; + dug=ug_epl_sub_iter->Duplicate();_assert_(dug); + ug_epl_sub_iter->Copy(dug);_assert_(ug_epl_sub_iter); + dug->AYPX(ug_epl,-1.0); + ndu_epl=dug->Norm(NORM_TWO); + delete dug; + nu_epl=ug_epl_sub_iter->Norm(NORM_TWO); + + if (xIsNan(ndu_epl) || xIsNan(nu_epl)) _error_("convergence criterion is NaN!"); + if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/ + if((ndu_epl/nu_epl)Duplicate();_assert_(ug_epl_sub_iter); + ug_epl->Copy(ug_epl_sub_iter);_assert_(ug_epl); + //end of the crap if(eplconverged){