Changeset 18719


Ignore:
Timestamp:
10/31/14 11:35:16 (10 years ago)
Author:
bdef
Message:

NEW:shfting hydrology lock to the nodes

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

Legend:

Unmodified
Added
Removed

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

    r18705 r18719  
    3030        bool   isefficientlayer;
    3131        int    hydrology_model;
     32
    3233
    3334        /*Now, do we really want DC?*/
    … …  
    6162                iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
    6263        }
    63 
    64 }/*}}}*/
     64}/*}}}*/
     65
    6566void HydrologyDCEfficientAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
    6667
    … …  
    9899}/*}}}*/
    99100
     101void HydrologyDCEfficientAnalysis::InitZigZagCounter(FemModel* femmodel){/*{{{*/
     102        int*   eplzigzag_counter =NULL;
     103       
     104        eplzigzag_counter=xNewZeroInit<int>(femmodel->nodes->Size());
     105        femmodel->parameters->AddObject(new IntVecParam(EplZigZagCounterEnum,eplzigzag_counter,femmodel->nodes->Size()));
     106        xDelete<int>(eplzigzag_counter);
     107}/*}}}*/
     108
     109void HydrologyDCEfficientAnalysis::ResetCounter(FemModel* femmodel){/*{{{*/
     110
     111        int*     eplzigzag_counter=NULL;
     112        Element* element=NULL;
     113
     114        femmodel->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
     115        for(int i=0;i<femmodel->nodes->Size();i++){
     116
     117                eplzigzag_counter[i]=0;
     118        }
     119        femmodel->parameters->SetParam(eplzigzag_counter,femmodel->nodes->Size(),EplZigZagCounterEnum);
     120        xDelete<int>(eplzigzag_counter);
     121}/*}}}*/
     122
    100123/*Finite Element Analysis*/
    101124void HydrologyDCEfficientAnalysis::Core(FemModel* femmodel){/*{{{*/
    102125        _error_("not implemented");
    103126}/*}}}*/
    104 ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/
     127
     128 ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/
    105129        /*Default, return NULL*/
    106130        return NULL;
    107131}/*}}}*/
    108 ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
     132
     133 ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
    109134_error_("Not implemented");
    110135}/*}}}*/
    111 ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/
     136
     137 ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/
    112138
    113139        /*Intermediaries*/
    … …  
    612638        xDelete<IssmDouble>(dbasis);
    613639}/*}}}*/
    614 void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask,Element* element){
     640void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, int* eplzigzag_counter, Element* element){
    615641
    616642        bool        active_element;
    … …  
    619645        IssmDouble  h_max;
    620646        IssmDouble  sedheadmin;
    621         Element*   basalelement=NULL;
     647        Element*    basalelement=NULL;
    622648
    623649        /*Get basal element*/
    … …  
    636662        /*Intermediaries*/
    637663
     664        int         penalty_lock;
    638665        int         numnodes      =basalelement->GetNumberOfNodes();
    639666        IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
    … …  
    642669        IssmDouble* eplhead       =xNew<IssmDouble>(numnodes);
    643670        IssmDouble* residual      =xNew<IssmDouble>(numnodes);
    644 
     671       
    645672        IssmDouble init_thick    =basalelement->GetMaterialParameter(HydrologydcEplInitialThicknessEnum);
    646673        IssmDouble colapse_thick =basalelement->GetMaterialParameter(HydrologydcEplColapseThicknessEnum);
    … …  
    648675        Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
    649676        active_element_input->GetInputValue(&active_element);
     677
     678        basalelement->parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum);
    650679
    651680        basalelement-> GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum);
    … …  
    662691                /*Activate EPL if residual is >0 */
    663692                if(residual[i]>0.){
     693                        if(old_active[i]==0) eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
    664694                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
    665695                }
    … …  
    670700                        /*If epl thickness gets under 10-3 initial thickness, close the layer*/
    671701                        if(epl_thickness[i]<colapse_thick){
    672                                 vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
    673                                 epl_thickness[i]=init_thick;
     702                                eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
     703                                if(eplzigzag_counter[basalelement->nodes[i]->Sid()]<penalty_lock |penalty_lock==0){
     704                                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
     705                                        epl_thickness[i]=init_thick;
     706                                }
     707                                else{
     708                                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
     709                                }
    674710                        }
    675711                }
    … …  
    678714
    679715                if(eplhead[i]>=h_max && active_element){
    680 
    681716                        for(j=0;j<numnodes;j++){
    682717                                if(old_active[j]>0.){
    … …  
    685720                                /*Increase of the domain is on the downstream node in term of sediment head*/
    686721                                if(sedhead[j] == sedheadmin){
     722                                        if(old_active[i]==0) eplzigzag_counter[basalelement->nodes[j]->Lid()] ++;
    687723                                        vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL);
    688724                                }
    689725                        }
    690726                }
    691                
    692         }
     727        }
     728        basalelement->AddInput(HydrologydcEplThicknessEnum,epl_thickness,basalelement->GetElementType());
     729
    693730        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    694731        xDelete<IssmDouble>(epl_thickness);
    … …  
    723760               
    724761        basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum);
     762
     763        bool active_element;
     764       
     765        Input*  active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);           
     766        active_element_input->GetInputValue(&active_element);
     767
    725768       
    726769        for(int i=0;i<numnodes;i++) flag+=active[i];
    … …  
    730773                        active_vec->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
    731774                }
     775        }
     776        else if(active_element){
     777                /*Checking Stuff*/
     778                for(int i=0;i<numnodes;i++){
     779                        active_vec->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
     780                }               
    732781        }
    733782        else{

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

    r18645 r18719  
    2020                void CreateConstraints(Constraints* constraints,IoModel* iomodel);
    2121                void CreateLoads(Loads* loads, IoModel* iomodel);
    22 
     22                void InitZigZagCounter(FemModel* femmodel);
     23                void ResetCounter(FemModel* femmodel);
     24                       
    2325                /*Finite element Analysis*/
    2426                void           Core(FemModel* femmodel);
    … …  
    3840                IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
    3941                IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
    40                 void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask,Element* element);
     42                void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, int* eplzigzag_counter, Element* element);
    4143                void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element);
    4244                void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode);

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

    r18713 r18719  
    1717        int         penalty_lock;
    1818        int         hydro_maxiter;
    19         int*        elementactive_counter =NULL;
     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         if(isefficientlayer==1){
    60                 elementactive_counter=xNewZeroInit<int>(iomodel->numberofelements);
    61                 parameters->AddObject(new IntVecParam(ElementActiveCounterEnum,elementactive_counter,iomodel->numberofelements));
    62                 xDelete<int>(elementactive_counter);
    63         }
    6458       
    6559}/*}}}*/
    … …  
    10498                iomodel->FetchDataToInput(elements,HydrologydcMaskEplactiveNodeEnum);
    10599
    106                 for(int i=0;i<elements->Size();i++){
    107                         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    108                         Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input);
     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); */
    109103                       
    110                         if(node_mask_input->Max()>0.) {
    111                                 element_active = true;
    112                         }
    113                         else{
    114                                 element_active = false;
    115                         }
    116                         element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active));
    117                 }
     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                /* } */
    118112        }
    119113}/*}}}*/
    … …  
    684678
    685679        bool     element_active;
    686         bool     old_active;
    687         int      penalty_lock;
    688         int*     elementactive_counter=NULL;
    689680        Element* element=NULL;
    690681
    691         femmodel->parameters->FindParam(&elementactive_counter,NULL,ElementActiveCounterEnum);
    692682        for(int i=0;i<femmodel->elements->Size();i++){
    693683                element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    694                 element->parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum);
    695                
    696                 Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input);
    697                 Input* old_active_input = element->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(old_active_input);
    698                 old_active_input->GetInputValue(&old_active);
     684                        Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input);
    699685               
    700686                if(node_mask_input->Max()>0.){
    … …  
    704690                        element_active = false;
    705691                }
    706                 if(old_active!=element_active){
    707                         elementactive_counter[i]=elementactive_counter[i]+1;
    708                 }
    709                 if(elementactive_counter[i]>penalty_lock){
    710                         element_active = true;
    711                 }
    712692                element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active));
    713693        }
    714         femmodel->parameters->SetParam(elementactive_counter,femmodel->elements->Size(),ElementActiveCounterEnum);
    715         xDelete<int>(elementactive_counter);
    716 }/*}}}*/
    717 
    718 void HydrologyDCInefficientAnalysis::ResetCounter(FemModel* femmodel){/*{{{*/
    719 
    720         int*     elementactive_counter=NULL;
    721         Element* element=NULL;
    722         int N;
    723 
    724         femmodel->parameters->FindParam(&elementactive_counter,&N,ElementActiveCounterEnum);
    725         for(int i=0;i<femmodel->elements->Size();i++){
    726                 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    727                
    728                 elementactive_counter[i]=0;
    729         }               
    730         femmodel->parameters->SetParam(elementactive_counter,femmodel->elements->Size(),ElementActiveCounterEnum);
    731         xDelete<int>(elementactive_counter);
    732 }/*}}}*/
     694}/*}}}*/

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

    r18705 r18719  
    4141                IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
    4242                void ElementizeEplMask(FemModel* femmodel);
    43                 void ResetCounter(FemModel* femmodel);
    4443};
    4544#endif

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

    r18715 r18719  
    18091809void FemModel::HydrologyEPLupdateDomainx(void){ /*{{{*/
    18101810
    1811         Vector<IssmDouble>* mask          = NULL;
    1812         IssmDouble*         serial_mask   = NULL;
    1813         Vector<IssmDouble>* active        = NULL;
    1814         IssmDouble*         serial_active = NULL;
    1815        
     1811        Vector<IssmDouble>* mask                                                        = NULL;
     1812        IssmDouble*         serial_mask                         = NULL;
     1813        Vector<IssmDouble>* active                                              = NULL;
     1814        IssmDouble*         serial_active                       = NULL;
     1815        int*                eplzigzag_counter = NULL;
     1816               
    18161817        HydrologyDCEfficientAnalysis* effanalysis =  new HydrologyDCEfficientAnalysis();
    18171818        /*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/
    1818         mask=new Vector<IssmDouble>(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
     1819        mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
     1820        this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
    18191821
    18201822        for (int i=0;i<elements->Size();i++){
    18211823                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    1822                 effanalysis->HydrologyEPLGetMask(mask,element);
    1823         }
    1824 
     1824                effanalysis->HydrologyEPLGetMask(mask,eplzigzag_counter,element);
     1825        }
     1826        printf("POINTER = %p\n",eplzigzag_counter);
     1827        this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum);
    18251828        /*Assemble and serialize*/
    18261829        mask->Assemble();
    18271830        serial_mask=mask->ToMPISerial();
     1831        xDelete<int>(eplzigzag_counter);
    18281832        delete mask;
    18291833

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r18717 r18719  
    109109        EplHeadSlopeXEnum,
    110110        EplHeadSlopeYEnum,
    111         ElementActiveCounterEnum,
     111        EplZigZagCounterEnum,
    112112        HydrologydcMaxIterEnum,
    113113        HydrologydcRelTolEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r18717 r18719  
    117117                case EplHeadSlopeXEnum : return "EplHeadSlopeX";
    118118                case EplHeadSlopeYEnum : return "EplHeadSlopeY";
    119                 case ElementActiveCounterEnum : return "ElementActiveCounter";
     119                case EplZigZagCounterEnum : return "EplZigZagCounter";
    120120                case HydrologydcMaxIterEnum : return "HydrologydcMaxIter";
    121121                case HydrologydcRelTolEnum : return "HydrologydcRelTol";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r18718 r18719  
    117117              else if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum;
    118118              else if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum;
    119               else if (strcmp(name,"ElementActiveCounter")==0) return ElementActiveCounterEnum;
     119              else if (strcmp(name,"EplZigZagCounter")==0) return EplZigZagCounterEnum;
    120120              else if (strcmp(name,"HydrologydcMaxIter")==0) return HydrologydcMaxIterEnum;
    121121              else if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum;

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

    r18705 r18719  
    6969                /*Initialize the element mask*/
    7070                inefanalysis->ElementizeEplMask(femmodel);
     71                effanalysis->InitZigZagCounter(femmodel);
    7172        }
    7273        /*The real computation starts here, outermost loop is on the two layer system*/
    … …  
    222223                                        InputUpdateFromConstantx(femmodel,eplconverged,ConvergedEnum);
    223224                                        InputUpdateFromSolutionx(femmodel,ug_epl);
    224                                         inefanalysis->ResetCounter(femmodel);
     225                                        effanalysis->ResetCounter(femmodel);
    225226                                        break;
    226227                                }
Note: See TracChangeset for help on using the changeset viewer.