Changeset 16954


Ignore:
Timestamp:
11/26/13 13:43:30 (11 years ago)
Author:
bdef
Message:

CHG: Changing the convergance mecanism in hydro

Location:
issm/trunk-jpl/src/c
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

    r16899 r16954  
    392392                int numvertices = element->GetNumberOfVertices();
    393393                Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum);
    394                 if(node_mask_input->Max()>0.) element_active = true;
    395                 else                          element_active = false;
    396 
     394                if(node_mask_input->Max()>0.) {
     395                        element_active = true;
     396                }
     397                else{
     398                        element_active = false;
     399                }
    397400                element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active));
    398401        }

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r16949 r16954  
    276276
    277277                #ifdef _HAVE_HYDROLOGY_
     278                virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index)=0;
    278279                virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0;
    279                 virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, int index)=0;
    280280                virtual void GetHydrologyTransfer(Vector<IssmDouble>* transfer)=0;
    281281                virtual void HydrologyEPLGetMask(Vector<IssmDouble>* mask)=0;

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r16947 r16954  
    65446544                transfer_input->GetInputValue(&transfer,gauss);
    65456545                scalar = Jdet*gauss->weight*(water_load+transfer);
    6546                 //printf("are we loading: load,%e ,transfer,%e\n",water_load, transfer);
    65476546                if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
    65486547                for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
    … …  
    67266725/*}}}*/
    67276726/*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/
    6728 void  Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){
    6729 
    6730                 int        hmax_flag;
    6731                 IssmDouble h_max;
    6732                 IssmDouble rho_ice,rho_water;
    6733                 IssmDouble thickness,bed;
    6734                 /*Get the flag to the limitation method*/
    6735                 this->parameters->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
    6736 
    6737                 /*Switch between the different cases*/
    6738                 switch(hmax_flag){
    6739                         case 0:
    6740                                 h_max=1.0e+10;
    6741                                 break;
    6742                         case 1:
    6743                                 parameters->FindParam(&h_max,HydrologydcSedimentlimitEnum);
    6744                                 break;
    6745                         case 2:
    6746                                 rho_ice=matpar->GetRhoIce();
    6747                                 rho_water=matpar->GetRhoFreshwater();
    6748                                 this->GetInputValue(&thickness,innode,ThicknessEnum);
    6749                                 this->GetInputValue(&bed,innode,BedEnum);
    6750                                 h_max=((rho_ice*thickness)/rho_water)+bed;
    6751                                 break;
    6752                         case 3:
    6753                                 _error_("Using normal stress  not supported yet");
    6754                                 break;
    6755                         default:
    6756                                 _error_("no case higher than 3 for SedimentlimitFlag");
    6757                 }
    6758                 /*Assign output pointer*/
    6759                 *ph_max=h_max;
    6760 }
    6761 /*}}}*/
    6762 /*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/
    67636727void  Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){
    6764 
     6728       
    67656729        int        hmax_flag;
    67666730        IssmDouble h_max;
    … …  
    67696733        /*Get the flag to the limitation method*/
    67706734        this->parameters->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
    6771 
     6735       
    67726736        /*Switch between the different cases*/
    67736737        switch(hmax_flag){
    6774                 case 0:
    6775                         h_max=1.0e+10;
    6776                         break;
    6777                 case 1:
    6778                         parameters->FindParam(&h_max,HydrologydcSedimentlimitEnum);
    6779                         break;
    6780                 case 2:
    6781                         rho_ice=matpar->GetRhoIce();
    6782                         rho_water=matpar->GetRhoFreshwater();
    6783                         this->GetInputValue(&thickness,this->nodes[index],ThicknessEnum);
    6784                         this->GetInputValue(&bed,this->nodes[index],BedEnum);
    6785                         h_max=((rho_ice*thickness)/rho_water)+bed;
    6786                         break;
    6787                 case 3:
    6788                         _error_("Using normal stress  not supported yet");
    6789                         break;
    6790                 default:
    6791                         _error_("no case higher than 3 for SedimentlimitFlag");
     6738        case 0:
     6739                h_max=1.0e+10;
     6740                break;
     6741        case 1:
     6742                parameters->FindParam(&h_max,HydrologydcSedimentlimitEnum);
     6743                break;
     6744        case 2:
     6745                rho_ice=matpar->GetRhoIce();
     6746                rho_water=matpar->GetRhoFreshwater();
     6747                this->GetInputValue(&thickness,this->nodes[index],ThicknessEnum);
     6748                this->GetInputValue(&bed,this->nodes[index],BedEnum);
     6749                h_max=((rho_ice*thickness)/rho_water)+bed;
     6750                break;
     6751        case 3:
     6752                _error_("Using normal stress  not supported yet");
     6753                break;
     6754        default:
     6755                _error_("no case higher than 3 for SedimentlimitFlag");
     6756        }
     6757        /*Assign output pointer*/
     6758        *ph_max=h_max;
     6759}
     6760/*}}}*/
     6761/*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/
     6762void  Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){
     6763       
     6764        int        hmax_flag;
     6765        IssmDouble h_max;
     6766        IssmDouble rho_ice,rho_water;
     6767        IssmDouble thickness,bed;
     6768        /*Get the flag to the limitation method*/
     6769        this->parameters->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
     6770       
     6771        /*Switch between the different cases*/
     6772        switch(hmax_flag){
     6773        case 0:
     6774                h_max=1.0e+10;
     6775                break;
     6776        case 1:
     6777                parameters->FindParam(&h_max,HydrologydcSedimentlimitEnum);
     6778                break;
     6779        case 2:
     6780                rho_ice=matpar->GetRhoIce();
     6781                rho_water=matpar->GetRhoFreshwater();
     6782                this->GetInputValue(&thickness,innode,ThicknessEnum);
     6783                this->GetInputValue(&bed,innode,BedEnum);
     6784                h_max=((rho_ice*thickness)/rho_water)+bed;
     6785                break;
     6786        case 3:
     6787                _error_("Using normal stress  not supported yet");
     6788                break;
     6789        default:
     6790                _error_("no case higher than 3 for SedimentlimitFlag");
    67926791        }
    67936792        /*Assign output pointer*/
    … …  
    68316830                        active_element_input->GetInputValue(&active_element);
    68326831
    6833                         GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveNodeEnum);
    68346832                        GetInputListOnVertices(&sed_head[0],SedimentHeadEnum);
    68356833                        GetInputListOnVertices(&epl_head[0],EplHeadEnum);
    … …  
    68616859                                               
    68626860                                                /*No transfer if the sediment head is allready at the maximum*/
    6863                                                 this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
     6861                                                this->GetHydrologyDCInefficientHmax(&h_max,i);
    68646862                                                if(sed_head[i]>=h_max)wh_trans=0.0;
    68656863                                        }
    … …  
    68706868                                        /*Assign output pointer*/
    68716869                                        transfer->SetValue(doflist[i],wh_trans,INS_VAL);
     6870                                        /* if(nodes[i]->id>=54){ */
     6871                                        /*      printf("%i %e %e %e \n",nodes[i]->id-54,wh_trans,sed_head[i],epl_head[i]); */
     6872                                        /* } */
     6873                                        /* else{*/
     6874                                        /*      printf("%i %e %e %e \n",nodes[i]->id,wh_trans,sed_head[i],epl_head[i]); */
    68726875                                }
    68736876                        }
    … …  
    68906893
    68916894        GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum);
    6892 
     6895       
    68936896        for(int i=0;i<numnodes;i++) flag+=active[i];
    68946897
    … …  
    69016904                /*Do not do anything: at least one node is active for this element but this element is not solved for*/
    69026905        }
    6903 
    69046906}
    69056907/*}}}*/
    … …  
    69396941                }
    69406942                /*If epl thickness gets under 0, close the layer*/
    6941                 else if(epl_thickness[i]<0.0){
    6942                         vec_mask->SetValue(nodes[i]->Sid(),0.,INS_VAL);
    6943                 }
     6943                /* else if(epl_thickness[i]<0.0){ */
     6944                /*      vec_mask->SetValue(nodes[i]->Sid(),0.,INS_VAL); */
     6945                /* } */
    69446946                /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
    6945                 this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
     6947                this->GetHydrologyDCInefficientHmax(&h_max,i);
    69466948                if(eplhead[i]>=h_max && this->AnyActive()){
    69476949                        for(j=0;j<numdof;j++){

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r16949 r16954  
    343343                ElementVector* CreatePVectorL2ProjectionEPL(void);
    344344                void           CreateHydrologyWaterVelocityInput(void);
     345                void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index);
    345346                void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    346                 void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index);
    347347                void           GetHydrologyTransfer(Vector<IssmDouble>* transfer);
    348348                void           HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);

  • issm/trunk-jpl/src/c/cores/hydrology_core.cpp

    r16887 r16954  
    9494                                if(VerboseSolution()) _printf0_("   saving results \n");
    9595                                if(isefficientlayer){
    96                                         int outputs[8] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveNodeEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum,HydrologydcEplThicknessEnum,WaterTransferEnum};
    97                                         femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],8);
     96                                        int outputs[9] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveNodeEnum,HydrologydcMaskEplactiveEltEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum,HydrologydcEplThicknessEnum,WaterTransferEnum};
     97                                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],9);
    9898                                }
    9999                                else{

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp

    r16918 r16954  
    2020        Vector<IssmDouble>* uf_epl=NULL;
    2121        Vector<IssmDouble>* uf_epl_sub_iter=NULL;
     22        Vector<IssmDouble>* ug_epl_sub_iter=NULL;
    2223        Vector<IssmDouble>* ug_epl_main_iter=NULL;
     24
    2325
    2426        Vector<IssmDouble>* ys=NULL;
    2527        Vector<IssmDouble>* dug=NULL;
     28
     29        //testing stuff
     30        Vector<IssmDouble>* duf=NULL;
    2631
    2732        Matrix<IssmDouble>* Kff=NULL;
    … …  
    5257
    5358        /*Retrieve inputs as the initial state for the non linear iteration*/
    54         GetSolutionFromInputsx(&ug_sed,femmodel);
     59        GetSolutionFromInputsx(&ug_sed,femmodel);       
     60
     61        //test
     62        GetSolutionFromInputsx(&uf_sed,femmodel);_assert_(uf_sed);
    5563
    5664        if(isefficientlayer) {
    5765                femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
    5866                GetSolutionFromInputsx(&ug_epl,femmodel);
     67
     68                //test
     69                GetSolutionFromInputsx(&uf_epl,femmodel);_assert_(uf_epl);
     70
    5971                /*Initialize the transfer input*/
    6072                HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis();
    … …  
    7284                ug_sed_main_iter=ug_sed->Duplicate();
    7385                ug_sed->Copy(ug_sed_main_iter);
     86               
     87                //test
     88                uf_sed_sub_iter=uf_sed->Duplicate();
     89                uf_sed->Copy(uf_sed_sub_iter);
     90
    7491                if(isefficientlayer){
    7592                        ug_epl_main_iter=ug_epl->Duplicate();
    7693                        ug_epl->Copy(ug_epl_main_iter);
     94                        //test
     95                        ug_epl_sub_iter=ug_epl->Duplicate();
     96                        ug_epl->Copy(ug_epl_sub_iter);
     97
    7798                }
    7899
    … …  
    95116                        Solverx(&uf_sed,Kff,pf,uf_sed_sub_iter,df,femmodel->parameters);
    96117                        delete Kff; delete pf; delete df;
    97                         delete uf_sed_sub_iter;
    98                         uf_sed_sub_iter=uf_sed->Duplicate();
    99                         uf_sed->Copy(uf_sed_sub_iter);
     118                        /* delete uf_sed_sub_iter; */
     119                        /* uf_sed_sub_iter=uf_sed->Duplicate(); */
     120                        /* uf_sed->Copy(uf_sed_sub_iter); */
    100121                        delete ug_sed;
    101122                        Mergesolutionfromftogx(&ug_sed,uf_sed,ys,femmodel->nodes,femmodel->parameters); delete ys;
    … …  
    110131                                }
    111132                        }
     133
    112134                        sedcount++;
     135
     136                        //testing stuff
     137                        if(sedconverged){
     138                                sedconverged=false;
     139                                duf=uf_sed_sub_iter->Duplicate();_assert_(duf);
     140                                uf_sed_sub_iter->Copy(duf);_assert_(uf_sed_sub_iter);
     141                                duf->AYPX(uf_sed,-1.0);
     142                                ndu_sed=duf->Norm(NORM_TWO);
     143                                delete duf;
     144                                nu_sed=uf_sed_sub_iter->Norm(NORM_TWO);
     145                                if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("convergence criterion is NaN!");
     146                                if((ndu_sed/nu_sed)<eps_hyd){
     147                                if(VerboseConvergence()) _printf0_("   # Inner sediment convergence achieve \n");
     148                                        sedconverged=true;
     149                                }
     150                        }
     151                        delete uf_sed_sub_iter;
     152                        uf_sed_sub_iter=uf_sed->Duplicate();_assert_(uf_sed_sub_iter);
     153                        uf_sed->Copy(uf_sed_sub_iter);_assert_(uf_sed);
     154                        //end of the crap
    113155
    114156                        if(sedconverged){
    … …  
    163205                                Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters);
    164206                                delete Kff; delete pf; delete df;
    165                                 delete uf_epl_sub_iter; 
     207                                delete uf_epl_sub_iter;
    166208                                uf_epl_sub_iter=uf_epl->Duplicate();
    167209                                uf_epl->Copy(uf_epl_sub_iter);
    … …  
    171213                                ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
    172214                                femmodel->HydrologyEPLupdateDomainx();                 
     215                                        /* /\*Updating Nodal Mask*\/ */
     216                                        /* HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); */
     217                                        /* analysis->ElementizeEplMask(femmodel); */
     218                                        /* delete analysis; */
     219                                        /* femmodel->HydrologyTransferx(); */
    173220
    174221                                if (!eplconverged){
    … …  
    182229                                }
    183230                                eplcount++;
     231
     232                                //testing stuff
     233                                if(eplconverged){
     234                                        eplconverged=false;
     235                                        dug=ug_epl_sub_iter->Duplicate();_assert_(dug);
     236                                        ug_epl_sub_iter->Copy(dug);_assert_(ug_epl_sub_iter);
     237                                        dug->AYPX(ug_epl,-1.0);
     238                                        ndu_epl=dug->Norm(NORM_TWO);
     239                                        delete dug;
     240                                        nu_epl=ug_epl_sub_iter->Norm(NORM_TWO);
     241                                       
     242                                        if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("convergence criterion is NaN!");
     243                                        if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/
     244                                        if((ndu_epl/nu_epl)<eps_hyd)eplconverged=true;
     245                                }
     246                                delete ug_epl_sub_iter;
     247                                ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter);
     248                                ug_epl->Copy(ug_epl_sub_iter);_assert_(ug_epl);
     249                                //end of the crap
    184250
    185251                                if(eplconverged){
Note: See TracChangeset for help on using the changeset viewer.