Changeset 18975


Ignore:
Timestamp:
01/05/15 11:21:22 (10 years ago)
Author:
bdef
Message:

minor clean-up in Hydro Code

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

Legend:

Unmodified
Added
Removed

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

    r18903 r18975  
    2424        iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    2525
     26        /*If not return*/
    2627        if(!isefficientlayer) return;
     28
     29        /*If yes, initialize a flip flop counter*/
    2730        iomodel->FetchData(&eplflip_lock,HydrologydcEplflipLockEnum);
    28  
    2931        parameters->AddObject(new IntParam(HydrologydcEplflipLockEnum,eplflip_lock));
    3032       
    31         /*Nothing for now*/
    3233}/*}}}*/
    3334void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    … …  
    3536        bool   isefficientlayer;
    3637        int    hydrology_model;
    37 
    3838
    3939        /*Now, do we really want DC?*/
    … …  
    6161        iomodel->FetchDataToInput(elements,HydrologydcEplInitialThicknessEnum);
    6262        iomodel->FetchDataToInput(elements,HydrologydcEplMaxThicknessEnum);
    63         iomodel->FetchDataToInput(elements,HydrologydcSedimentTransmitivityEnum);
    6463        iomodel->FetchDataToInput(elements,HydrologydcEplThicknessEnum);
    65         if(iomodel->domaintype!=Domain2DhorizontalEnum){
     64                if(iomodel->domaintype!=Domain2DhorizontalEnum){
    6665                iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
    6766                iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
    … …  
    8584        iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
    8685}/*}}}*/
     86
    8787void HydrologyDCEfficientAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
    8888
    … …  
    9898
    9999        IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum,P1Enum);
    100 
    101 }/*}}}*/
     100}/*}}}*/
     101
    102102void HydrologyDCEfficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
    103103        /*Nothing for now*/
    … …  
    131131}/*}}}*/
    132132
    133  ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/
     133ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/
    134134        /*Default, return NULL*/
    135135        return NULL;
    136136}/*}}}*/
    137137
    138  ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
     138ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
    139139_error_("Not implemented");
    140140}/*}}}*/
    141141
    142  ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/
     142ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/
    143143
    144144        /*Intermediaries*/
    … …  
    244244        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    245245        return Ke;
    246 
    247 }/*}}}*/
     246}/*}}}*/
     247
    248248ElementVector* HydrologyDCEfficientAnalysis::CreatePVector(Element* element){/*{{{*/
    249249
    … …  
    348348        return pe;
    349349}/*}}}*/
     350
    350351void HydrologyDCEfficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    351352        element->GetSolutionFromInputsOneDof(solution,EplHeadEnum);
    352353}/*}}}*/
     354
    353355void HydrologyDCEfficientAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
    354356        _error_("Not implemented yet");
    355357}/*}}}*/
     358
    356359void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    357360
    … …  
    375378        /*Fetch number of nodes for this finite element*/
    376379        int numnodes = basalelement->GetNumberOfNodes();
    377 
    378380
    379381        /*Fetch dof list and allocate solution vector*/
    … …  
    401403        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    402404} /*}}}*/
     405
    403406void HydrologyDCEfficientAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
    404407        /*Default, do nothing*/
    … …  
    415418        return rho_freshwater*g*epl_porosity*(water_compressibility+(epl_compressibility/epl_porosity));                 
    416419}/*}}}*/
     420
    417421IssmDouble HydrologyDCEfficientAnalysis::SedimentStoring(Element* element){/*{{{*/
    418422        IssmDouble rho_freshwater           = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    … …  
    535539                _error_("not Implemented Yet");
    536540                }
    537                        
     541               
    538542                int         numnodes      = element->GetNumberOfNodes();
    539543                IssmDouble* thickness     = xNew<IssmDouble>(numnodes);
    … …  
    602606                xDelete<IssmDouble>(bed);
    603607        }
    604 }
    605 /*}}}*/
     608}/*}}}*/
     609
    606610void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    607611        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    … …  
    631635        xDelete<IssmDouble>(dbasis);
    632636}/*}}}*/
     637
    633638void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, int* eplzigzag_counter, Element* element){
    634639
    … …  
    655660        /*Intermediaries*/
    656661
    657         int         penalty_lock;
     662        int         eplflip_lock;
    658663        int         numnodes      =basalelement->GetNumberOfNodes();
    659664        IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
    … …  
    669674        active_element_input->GetInputValue(&active_element);
    670675
    671         basalelement->parameters->FindParam(&penalty_lock,HydrologydcEplflipLockEnum);
     676        basalelement->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum);
    672677
    673678        basalelement-> GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum);
    … …  
    691696                else if(old_active[i]>0.){
    692697                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
    693                         /*If epl thickness gets under 10-3 initial thickness, close the layer*/
     698                        /*If epl thickness gets under colapse thickness, close the layer*/
    694699                        if(epl_thickness[i]<colapse_thick){
    695700                                eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
    696                                 if(eplzigzag_counter[basalelement->nodes[i]->Lid()]<penalty_lock |penalty_lock==0){
     701                                /*Avoid flipfloping between open and closed states*/
     702                                if(eplzigzag_counter[basalelement->nodes[i]->Lid()]<eplflip_lock |eplflip_lock==0){
    697703                                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
    698704                                        epl_thickness[i]=init_thick;
    … …  
    705711                /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
    706712                GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->nodes[i]);
    707 
    708713                if(eplhead[i]>=h_max && active_element){
    709714                        for(j=0;j<numnodes;j++){
    … …  
    729734}
    730735/*}}}*/
     736
    731737void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){
    732738        /*Constants*/
    … …  
    750756        IssmDouble  flag     = 0.;
    751757        IssmDouble* active   = xNew<IssmDouble>(numnodes);
    752                
     758
     759        /*Pass the activity mask from elements to nodes*/
    753760        basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum);
    754761        bool active_element;
    … …  
    758765        for(int i=0;i<numnodes;i++) flag+=active[i];
    759766
     767        /*If any node is active all the node in the element are active*/
    760768        if(flag>0.){
    761769                for(int i=0;i<numnodes;i++){
    … …  
    763771                }
    764772        }
     773        /*If the element is active all its nodes are active*/
    765774        else if(active_element){
    766775                for(int i=0;i<numnodes;i++){

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

    r18903 r18975  
    1010        return 1;
    1111}/*}}}*/
     12
    1213void HydrologyDCInefficientAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    1314
    … …  
    1718        int         penalty_lock;
    1819        int         hydro_maxiter;
    19         /* int*        elementactive_counter =NULL; */
    2020        bool        isefficientlayer;
    2121        IssmDouble  sedimentlimit;
    … …  
    5656        parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock));
    5757        parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter));
    58        
    59 }/*}}}*/
     58}/*}}}*/
     59
    6060void HydrologyDCInefficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    6161
    … …  
    9797        if(isefficientlayer){
    9898                iomodel->FetchDataToInput(elements,HydrologydcMaskEplactiveNodeEnum);
    99 
    100                 /* for(int i=0;i<elements->Size();i++){ */
    101                 /*      Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i)); */
    102                 /*      Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input); */
    103                        
    104                 /*      if(node_mask_input->Max()>0.) { */
    105                 /*              element_active = true; */
    106                 /*      } */
    107                 /*      else{ */
    108                 /*              element_active = false; */
    109                 /*      } */
    110                 /*      element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active)); */
    111                 /* } */
    112         }
    113 }/*}}}*/
     99        }
     100}/*}}}*/
     101
    114102void HydrologyDCInefficientAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
    115103
    … …  
    125113        iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
    126114}/*}}}*/
     115
    127116void HydrologyDCInefficientAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
    128117
    … …  
    134123        IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyDCInefficientAnalysisEnum,P1Enum);
    135124}/*}}}*/
     125
    136126void HydrologyDCInefficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
    137127
    … …  
    165155        _error_("not implemented");
    166156}/*}}}*/
     157
    167158ElementVector* HydrologyDCInefficientAnalysis::CreateDVector(Element* element){/*{{{*/
    168159        /*Default, return NULL*/
    169160        return NULL;
    170161}/*}}}*/
     162
    171163ElementMatrix* HydrologyDCInefficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
    172164_error_("Not implemented");
    173165}/*}}}*/
     166
    174167ElementMatrix* HydrologyDCInefficientAnalysis::CreateKMatrix(Element* element){/*{{{*/
    175168
    … …  
    275268        return Ke;
    276269}/*}}}*/
     270
    277271ElementVector* HydrologyDCInefficientAnalysis::CreatePVector(Element* element){/*{{{*/
    278272
    … …  
    312306        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    313307
    314 
    315308        /*Retrieve all inputs and parameters*/
    316309        basalelement->GetVerticesCoordinates(&xyz_list);
    … …  
    342335                /*Loading term*/
    343336                water_input->GetInputValue(&water_load,gauss);
    344 
    345337                scalar = Jdet*gauss->weight*(water_load);
    346 
    347338                if(dt!=0.) scalar = scalar*dt;
    348339                for(int i=0;i<numnodes;i++){
    … …  
    378369        return pe;
    379370}/*}}}*/
     371
    380372void HydrologyDCInefficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    381373        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    … …  
    405397        xDelete<IssmDouble>(dbasis);
    406398}/*}}}*/
     399
    407400void HydrologyDCInefficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    408401        element->GetSolutionFromInputsOneDof(solution,SedimentHeadEnum);
    409402}/*}}}*/
     403
    410404void HydrologyDCInefficientAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
    411405        _error_("Not implemented yet");
    412406}/*}}}*/
     407
    413408void HydrologyDCInefficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    414409
    … …  
    489484        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    490485}/*}}}*/
     486
    491487void HydrologyDCInefficientAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
    492488        /*Default, do nothing*/
    … …  
    551547        return h_max;
    552548}/*}}}*/
     549
    553550void  HydrologyDCInefficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/
    554551       

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r18953 r18975  
    19791979        counter=sum_counter;
    19801980        if(VerboseSolution()) _printf0_("   Number of active nodes L2 Projection: "<< counter <<"\n");
    1981 
    1982         /*+++++++this is done in the solution sequence+++++++++++++++*/
    1983         /*Update dof indexings*/
    1984         //      this->UpdateConstraintsx();
    1985         /*++++++++++++++++++++++++++*/
    1986 }
    1987 /*}}}*/
     1981}
     1982/*}}}*/

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

    r18726 r18975  
    127127               
    128128                        /*}}}*//*End of the sediment penalization loop*/
    129                         /*Update EPL mask*/
    130                        
    131                         /*++++++This is probably useless+++++++*/
    132                         /* if(isefficientlayer){ */
    133                         /*      inefanalysis->ElementizeEplMask(femmodel); */
    134                         /* } */
    135                         /*+++++++++++++*/
    136 
    137129                        sedconverged=false;
    138130                       
Note: See TracChangeset for help on using the changeset viewer.