Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 17161) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 17162) @@ -117,7 +117,5 @@ /*Intermediaries*/ int meshtype; - bool active_element; Element* basalelement; - Input* active_element_input=NULL; /*Get basal element*/ @@ -134,10 +132,6 @@ } - /* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */ - /* active_element_input->GetInputValue(&active_element); */ - /*Check that all nodes are active, else return empty matrix*/ if(!basalelement->AllActive()) { - //if(!active_element){ if(meshtype!=Mesh2DhorizontalEnum){ basalelement->DeleteMaterials(); @@ -210,8 +204,5 @@ /*Intermediaries*/ int meshtype; - bool active_element; Element* basalelement; - Input* active_element_input=NULL; - /*Get basal element*/ @@ -227,10 +218,6 @@ default: _error_("mesh "<GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */ - /* active_element_input->GetInputValue(&active_element); */ - /*Check that all nodes are active, else return empty matrix*/ if(!basalelement->AllActive()) { - //if(!active_element){ if(meshtype!=Mesh2DhorizontalEnum){ basalelement->DeleteMaterials(); @@ -312,16 +299,53 @@ void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - int meshtype; + + int meshtype,i; + Element* basalelement=NULL; + element->FindParam(&meshtype,MeshTypeEnum); - switch(meshtype){ - case Mesh2DhorizontalEnum: - element->InputUpdateFromSolutionOneDof(solution,EplHeadEnum); - break; - case Mesh3DEnum: - element->InputUpdateFromSolutionOneDofCollapsed(solution,EplHeadEnum); - break; - default: _error_("mesh "<IsOnBed()) return; + basalelement=element->SpawnBasalElement(); + } + else{ + basalelement = element; + } + + /*Intermediary*/ + int* doflist = NULL; + + /*Fetch number of nodes for this finite element*/ + int numnodes = basalelement->GetNumberOfNodes(); + + /*Fetch dof list and allocate solution vector*/ + basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + IssmDouble* eplHeads = xNew(numnodes); + IssmDouble* eplOldHeads = xNew(numnodes); + IssmDouble Stepping; + + + /*Use the dof list to index into the solution vector: */ + for(i=0;i(eplHeads[i])) _error_("NaN found in solution vector"); + } + + /*Get previous water head*/ + basalelement->GetInputListOnNodes(&eplOldHeads[0],EplHeadOldEnum); + + for(i=0;iAddBasalInput(EplHeadEnum,eplHeads,P1Enum); + + /*Free ressources:*/ + xDelete(eplHeads); + xDelete(eplOldHeads); + xDelete(doflist); + if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; +} /*}}}*/ /*Intermediaries*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17161) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17162) @@ -1566,6 +1566,6 @@ name==SedimentHeadOldEnum || name==SedimentHeadEnum || + name==EplHeadEnum || name==EplHeadOldEnum || - name==EplHeadEnum || name==HydrologydcEplThicknessOldEnum || name==HydrologydcEplInitialThicknessEnum || @@ -4725,16 +4725,18 @@ const int numdof = NDOF1 *NUMVERTICES; int *doflist = NULL; + int analysis_type; bool isefficientlayer; bool active_element; - int transfermethod,step; + int transfermethod; IssmDouble leakage,h_max; IssmDouble wh_trans,sed_thick; + IssmDouble epl_specificstoring,sedstoring; IssmDouble activeEpl[numdof],epl_thickness[numdof]; - IssmDouble epl_specificstoring[numdof],sedstoring[numdof]; IssmDouble epl_head[numdof],sed_head[numdof]; - IssmDouble old_transfer[numdof],sed_trans[numdof]; + IssmDouble preceding_transfer[numdof],sed_trans[numdof]; Input* active_element_input=NULL; + parameters->FindParam(&analysis_type,AnalysisTypeEnum); GetDofList(&doflist,NoneApproximationEnum,GsetEnum); @@ -4742,11 +4744,9 @@ /*Get the flag to know if the efficient layer is present*/ this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); - - this->parameters->FindParam(&step,StepEnum); if(isefficientlayer){ /*Also get the flag to the transfer method*/ this->parameters->FindParam(&transfermethod,HydrologydcTransferFlagEnum); - + /*Switch between the different transfer methods cases*/ switch(transfermethod){ @@ -4755,10 +4755,9 @@ break; case 1: - active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); active_element_input->GetInputValue(&active_element); - GetInputListOnVertices(&sed_head[0],SedimentHeadEnum); + GetInputListOnVertices(&sed_head[0],SedimentHeadEnum); GetInputListOnVertices(&sed_trans[0],HydrologydcSedimentTransmitivityEnum); GetInputListOnVertices(&epl_head[0],EplHeadEnum); @@ -4773,17 +4772,14 @@ } else{ + GetInputListOnVertices(&preceding_transfer[0],WaterTransferEnum); + sedstoring=matpar->GetSedimentStoring(); + epl_specificstoring=matpar->GetEplSpecificStoring(); + + for(int i=0;iGetHydrologyDCInefficientHmax(&h_max,i); - GetInputListOnVertices(&old_transfer[0],WaterTransferEnum); - - for(int i=0;iGetEplSpecificStoring(); - sedstoring[i]=matpar->GetSedimentStoring(); - 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[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); + wh_trans=epl_specificstoring*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){ @@ -4791,20 +4787,14 @@ } } - /*EPL head lower than sediment head, transfer from the sediment to the epl*/ + /* 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[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); + wh_trans=sedstoring*sed_trans[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick); } - - /*Introduce relaxation*/ - // wh_trans=old_transfer[i]+0.66*(wh_trans-old_transfer[i]); - - /* if(this->Id()==27){ */ - /* printf("old %e, transfer is %e ,eplhead %e, sedhead %e, factor %e \n",old_transfer[i],wh_trans,epl_head[i],sed_head[i],epl_specificstoring[i]*epl_thickness[i]*sed_trans[i]); */ - /* } */ - + + /*Relaxation stuff*/ + wh_trans=preceding_transfer[i]+0.8*(wh_trans-preceding_transfer[i]); + /*Assign output pointer*/ transfer->SetValue(doflist[i],wh_trans,INS_VAL); - - } } @@ -4846,4 +4836,5 @@ int i,j; const int numdof = NDOF1 *NUMVERTICES; + IssmDouble init_thick; IssmDouble h_max; IssmDouble sedheadmin; @@ -4857,4 +4848,6 @@ Input* active_element_input=NULL; + init_thick = matpar->GetEplInitialThickness(); + active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); active_element_input->GetInputValue(&active_element); @@ -4876,15 +4869,14 @@ } - /*If mask was alread one, keep one*/ + /*If mask was already one, keep one*/ else if(old_active[i]>0.){ vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL); - } - /*If epl thickness gets under 0, close the layer*/ - else if(epl_thickness[i]<0.0){ - vec_mask->SetValue(nodes[i]->Sid(),0.,INS_VAL); + /*If epl thickness gets under , close the layer*/ + if(epl_thickness[i]<0.001*init_thick){ + 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,i); - //if(eplhead[i]>=h_max && this->AnyActive()){ if(eplhead[i]>=h_max && active_element){ for(j=0;jSetValue(nodes[j]->Sid(),1.,INS_VAL); - // break; } } @@ -4908,4 +4899,5 @@ const int numdof = NDOF1 *NUMVERTICES; bool isefficientlayer; + bool active_element; IssmDouble n,A,dt,init_thick; IssmDouble rho_water,rho_ice; @@ -4913,9 +4905,9 @@ IssmDouble EPL_N,epl_conductivity; IssmDouble activeEpl[numdof],thickness[numdof]; - IssmDouble eplhead[numdof], old_thickness[numdof]; + IssmDouble eplhead[numdof],old_eplhead[numdof]; IssmDouble epl_slopeX[numdof],epl_slopeY[numdof]; + IssmDouble preceding_thickness[numdof],old_thickness[numdof]; IssmDouble ice_thickness[numdof],bed[numdof]; - bool active_element; Input* active_element_input=NULL; @@ -4925,4 +4917,5 @@ if(isefficientlayer){ + /*For now, assuming just one way to compute EPL thickness*/ rho_water = matpar->GetRhoWater(); @@ -4935,6 +4928,4 @@ A = material->GetAbar(); - // GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveNodeEnum); - active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); active_element_input->GetInputValue(&active_element); @@ -4947,25 +4938,25 @@ GetInputListOnVertices(&bed[0],BedEnum); - //if(!this->AnyActive()){ if(!active_element){ /*Keeping thickness to initial value if EPL is not active*/ - for(int i=0;iId()==27){ */ - /* printf("old %e, thickness is %e ,eplhead %e, slopeX %e, slopeY %e\n",old_thickness[i],thickness[i],eplhead[i],epl_slopeX[i],epl_slopeY[i]); */ - /*}*/ + /*Relaxation stuff*/ + thickness[i] = preceding_thickness[i]+0.8*(thickness[i]-preceding_thickness[i]); } } Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 17161) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 17162) @@ -1423,5 +1423,5 @@ void FemModel::HydrologyTransferx(void){ /*{{{*/ - Vector* transferg=NULL; + Vector* transferg=NULL; /*Vector allocation*/ @@ -1429,4 +1429,5 @@ for (int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); element->GetHydrologyTransfer(transferg); @@ -1440,4 +1441,5 @@ } /*}}}*/ + void FemModel::HydrologyEPLThicknessx(void){ /*{{{*/ Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17161) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17162) @@ -34,4 +34,5 @@ bool sedconverged,eplconverged,hydroconverged; + bool transfered; bool isefficientlayer; int constraints_converged; @@ -43,4 +44,5 @@ IssmDouble ndu_sed,nu_sed; IssmDouble ndu_epl,nu_epl; + IssmDouble SedConv,EplConv; /*Recover parameters: */ @@ -51,5 +53,5 @@ /*FIXME, hardcoded, put on an enum*/ - hydro_maxiter=150; + hydro_maxiter=100; hydrocount=1; @@ -78,4 +80,6 @@ sedcount=1; eplcount=1; + SedConv=1.0; + EplConv=1.0; /*If there is two layers we need an outer loop value to compute convergence*/ @@ -129,8 +133,14 @@ /*Update Elemental Mask and compute transfer*/ if(isefficientlayer){ - HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); - analysis->ElementizeEplMask(femmodel); - delete analysis; - femmodel->HydrologyTransferx(); + if(SedConv<0.3){ + HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); + analysis->ElementizeEplMask(femmodel); + delete analysis; + femmodel->HydrologyTransferx(); + transfered=true; + } + else{ + transfered=false; + } } /*Checking convegence on the value of the sediment head*/ @@ -144,4 +154,5 @@ if (ndu_sed==0.0 && nu_sed==0.0) nu_sed=1.0e-6; /*Hacking the case where the layer is empty*/ if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << " aiming lower than " << eps_hyd*100 << " %\n"); + SedConv = ndu_sed/nu_sed*100; if((ndu_sed/nu_sed)parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum); InputUpdateFromConstantx(femmodel,sedconverged,ConvergedEnum); @@ -161,5 +172,4 @@ /* {{{ *//*Now dealing with the EPL in the same way*/ if(isefficientlayer){ - femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); /*updating mask*/ @@ -170,10 +180,37 @@ eplconverged=false; + for(;;){ + ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter); ug_epl->Copy(ug_epl_sub_iter); + + if(EplConv<0.3){ + /* {{{ *//*Retriev the EPL head slopes and compute EPL Thickness*/ + if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); + femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum); + femmodel->UpdateConstraintsL2ProjectionEPLx(); + femmodel->parameters->SetParam(EplHeadSlopeXEnum,InputToL2ProjectEnum); + solutionsequence_linear(femmodel); + femmodel->parameters->SetParam(EplHeadSlopeYEnum,InputToL2ProjectEnum); + solutionsequence_linear(femmodel); + femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); + femmodel->HydrologyEPLThicknessx(); + + //updating mask after the computation of the epl thickness (Allow to close too thin EPL) + femmodel->HydrologyEPLupdateDomainx(); + /* }}} */ - /*Loop on the epl layer to deal with the penalization, have been removed TO CHECK*/ - + /*Update Elemental Mask and compute transfer*/ + HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); + analysis->ElementizeEplMask(femmodel); + delete analysis; + femmodel->HydrologyTransferx(); + transfered=true; + } + else{ + transfered=false; + } + SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); @@ -189,26 +226,5 @@ InputUpdateFromSolutionx(femmodel,ug_epl); ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); - - /* {{{ *//*Start by retrieving the EPL head slopes and compute EPL Thickness*/ - if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); - femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum); - femmodel->UpdateConstraintsL2ProjectionEPLx(); - femmodel->parameters->SetParam(EplHeadSlopeXEnum,InputToL2ProjectEnum); - solutionsequence_linear(femmodel); - femmodel->parameters->SetParam(EplHeadSlopeYEnum,InputToL2ProjectEnum); - solutionsequence_linear(femmodel); - femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); - femmodel->HydrologyEPLThicknessx(); - - //updating mask after the computation of the epl thickness (Allow to close too thin EPL) - femmodel->HydrologyEPLupdateDomainx(); - /* }}} */ - - /*Update Elemental Mask and compute transfer*/ - HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); - analysis->ElementizeEplMask(femmodel); - delete analysis; - femmodel->HydrologyTransferx(); - + dug=ug_epl_sub_iter->Duplicate();_assert_(dug); ug_epl_sub_iter->Copy(dug); @@ -220,12 +236,15 @@ 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(VerboseConvergence()) _printf0_(setw(50) << left << " Inner EPL Convergence criterion:" << ndu_epl/nu_epl*100 << " aiming lower than " << eps_hyd*100 << " %\n"); + EplConv=ndu_epl/nu_epl*100; if((ndu_epl/nu_epl)=hydro_maxiter){ _error_(" maximum number of EPL iterations (" << hydro_maxiter << ") exceeded"); } + eplcount++; delete ug_epl_sub_iter; - if(eplconverged){ + if(eplconverged && transfered){ + if(VerboseSolution()) _printf0_("eplconverged...\n"); InputUpdateFromConstantx(femmodel,eplconverged,ConvergedEnum); InputUpdateFromConstantx(femmodel,sediment_kmax,MeltingOffsetEnum); @@ -259,5 +278,5 @@ 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 (!xIsNan(eps_hyd)){ - if ((ndu_epl/nu_epl)