Changeset 17372


Ignore:
Timestamp:
03/03/14 11:36:17 (11 years ago)
Author:
bdef
Message:

CHG: moving hydrology related routines from the element to the analysis

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

Legend:

Unmodified
Added
Removed

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

    r17362 r17372  
    367367        return rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));               
    368368}/*}}}*/
    369 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss){/*{{{*/
    370         int        hmax_flag;
    371         IssmDouble h_max;
    372         IssmDouble rho_ice,rho_water;
    373         IssmDouble thickness,bed;
    374 
    375         /*Get the flag to the limitation method*/
    376         element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
    377        
    378         /*Switch between the different cases*/
    379         switch(hmax_flag){
    380         case 0:
    381                 h_max=1.0e+10;
    382                 break;
    383         case 1:
    384                 element->FindParam(&h_max,HydrologydcSedimentlimitEnum);
    385                 break;
    386         case 2:
    387                 rho_water= element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    388                 rho_ice= element->GetMaterialParameter(MaterialsRhoIceEnum);
    389                 element->GetInputValue(&thickness,gauss,ThicknessEnum);
    390                 element->GetInputValue(&bed,gauss,BedEnum);
    391                 h_max=((rho_ice*thickness)/rho_water)+bed;
    392                 break;
    393         case 3:
    394                 _error_("Using normal stress  not supported yet");
    395                 break;
    396         default:
    397                 _error_("no case higher than 3 for SedimentlimitFlag");
    398         }
    399         return h_max;
    400 }/*}}}*/
    401369IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss){/*{{{*/
    402370       
    … …  
    422390                /* get input */
    423391
    424                         element->GetInputValue(&epl_thickness,gauss,HydrologydcEplThicknessEnum);
     392                element->GetInputValue(&epl_thickness,gauss,HydrologydcEplThicknessEnum);
    425393                element->GetInputValue(&sed_head,gauss,SedimentHeadEnum);
    426394                element->GetInputValue(&epl_head,gauss,EplHeadEnum);
    … …  
    519487                }
    520488                       
    521                 int         numnodes = element->GetNumberOfNodes();
    522                 IssmDouble  thickness[numnodes];
    523                 IssmDouble  eplhead[numnodes];
    524                 IssmDouble  epl_slopeX[numnodes],epl_slopeY[numnodes];
    525                 IssmDouble  old_thickness[numnodes];
    526                 IssmDouble  ice_thickness[numnodes],bed[numnodes];
     489                int         numnodes      = element->GetNumberOfNodes();
     490                IssmDouble* thickness     = xNew<IssmDouble>(numnodes);
     491                IssmDouble* eplhead       = xNew<IssmDouble>(numnodes);
     492                IssmDouble* epl_slopeX    = xNew<IssmDouble>(numnodes);
     493                IssmDouble* epl_slopeY    = xNew<IssmDouble>(numnodes);
     494                IssmDouble* old_thickness = xNew<IssmDouble>(numnodes);
     495                IssmDouble* ice_thickness = xNew<IssmDouble>(numnodes);
     496                IssmDouble* bed           = xNew<IssmDouble>(numnodes);
    527497
    528498                Input*  active_element_input=element->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);               
    … …  
    574544                }
    575545                element->AddInput(HydrologydcEplThicknessEnum,thickness,P1Enum);
     546                xDelete<IssmDouble>(thickness);
     547                xDelete<IssmDouble>(eplhead);
     548                xDelete<IssmDouble>(epl_slopeX);
     549                xDelete<IssmDouble>(epl_slopeY);
     550                xDelete<IssmDouble>(old_thickness);
     551                xDelete<IssmDouble>(ice_thickness);
     552                xDelete<IssmDouble>(bed);
    576553        }
    577554}
    … …  
    604581        xDelete<IssmDouble>(dbasis);
    605582}/*}}}*/
     583void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask,Element* element){
     584
     585        bool        active_element;
     586        int         i,j;
     587        int         meshtype;
     588        IssmDouble  h_max;
     589        IssmDouble  sedheadmin;
     590        Element*   basalelement=NULL;
     591
     592        /*Get basal element*/
     593        element->FindParam(&meshtype,MeshTypeEnum);
     594        switch(meshtype){
     595                case Mesh2DhorizontalEnum:
     596                        basalelement = element;
     597                        break;
     598                case Mesh3DEnum:
     599                        if(!element->IsOnBed()) return;
     600                        basalelement = element->SpawnBasalElement();
     601                        break;
     602                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     603        }
     604
     605        /*Intermediaries*/
     606
     607        int         numnodes      =basalelement->GetNumberOfNodes();
     608        IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
     609        IssmDouble* old_active    =xNew<IssmDouble>(numnodes);
     610        IssmDouble* sedhead       =xNew<IssmDouble>(numnodes);
     611        IssmDouble* eplhead       =xNew<IssmDouble>(numnodes);
     612        IssmDouble* residual      =xNew<IssmDouble>(numnodes);
     613
     614        IssmDouble init_thick = basalelement->GetMaterialParameter(HydrologydcEplInitialThicknessEnum);
     615
     616        Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
     617        active_element_input->GetInputValue(&active_element);
     618
     619        basalelement-> GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum);
     620        basalelement->  GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum); 
     621        basalelement->  GetInputListOnVertices(&sedhead[0],SedimentHeadEnum);
     622        basalelement->  GetInputListOnVertices(&eplhead[0],EplHeadEnum);
     623        basalelement->  GetInputListOnVertices(&residual[0],SedimentHeadResidualEnum);
     624
     625        /*Get minimum sediment head of the element*/
     626        sedheadmin=sedhead[0];
     627        for(i=1;i<numnodes;i++) if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i];
     628
     629        for(i=0;i<numnodes;i++){
     630                /*Activate EPL if residual is >0 */
     631                if(residual[i]>0.){
     632                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
     633                }
     634
     635                /*If mask was already one, keep one*/
     636                else if(old_active[i]>0.){
     637                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
     638                        /*If epl thickness gets under , close the layer*/
     639                        if(epl_thickness[i]<0.001*init_thick){
     640                                vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
     641                                epl_thickness[i]=init_thick;
     642                        }
     643                }
     644                /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
     645                GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->nodes[i]);
     646                if(eplhead[i]>=h_max && active_element){
     647                        for(j=0;j<numnodes;j++){
     648                                if(old_active[j]>0.){
     649                                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
     650                                }
     651                                /*Increase of the domain is on the downstream node in term of sediment head*/
     652                                if(sedhead[j] == sedheadmin){
     653                                        vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL);
     654                                }
     655                        }
     656                }
     657        }
     658        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     659        xDelete<IssmDouble>(epl_thickness);
     660        xDelete<IssmDouble>(old_active);
     661        xDelete<IssmDouble>(sedhead);
     662        xDelete<IssmDouble>(eplhead);
     663        xDelete<IssmDouble>(residual);
     664}
     665/*}}}*/
     666void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){
     667        /*Constants*/
     668
     669        int      meshtype;
     670        Element*   basalelement=NULL;
     671
     672        /*Get basal element*/
     673        element->FindParam(&meshtype,MeshTypeEnum);
     674        switch(meshtype){
     675                case Mesh2DhorizontalEnum:
     676                        basalelement = element;
     677                        break;
     678                case Mesh3DEnum:
     679                        if(!element->IsOnBed()) return;
     680                        basalelement = element->SpawnBasalElement();
     681                        break;
     682                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     683        }
     684       
     685        const int   numnodes = basalelement->GetNumberOfNodes();
     686        IssmDouble  flag     = 0.;
     687        IssmDouble* active   = xNew<IssmDouble>(numnodes);
     688               
     689        basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum);
     690       
     691        for(int i=0;i<numnodes;i++) flag+=active[i];
     692
     693        if(flag>0.){
     694                for(int i=0;i<numnodes;i++){
     695                        active_vec->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
     696                }
     697        }
     698        else{
     699                /*Do not do anything: at least one node is active for this element but this element is not solved for*/
     700        }
     701        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     702        xDelete<IssmDouble>(active);
     703}
     704
     705void  HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/
     706       
     707        int        hmax_flag;
     708        IssmDouble h_max;
     709        IssmDouble rho_ice,rho_water;
     710        IssmDouble thickness,bed;
     711        /*Get the flag to the limitation method*/
     712        element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
     713       
     714        /*Switch between the different cases*/
     715        switch(hmax_flag){
     716        case 0:
     717                h_max=1.0e+10;
     718                break;
     719        case 1:
     720                element->FindParam(&h_max,HydrologydcSedimentlimitEnum);
     721                break;
     722        case 2:
     723                /*Compute max*/
     724                rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
     725                rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
     726                element->GetInputValue(&thickness,innode,ThicknessEnum);
     727                element->GetInputValue(&bed,innode,BedEnum);
     728                h_max=((rho_ice*thickness)/rho_water)+bed;
     729                break;
     730        case 3:
     731                _error_("Using normal stress  not supported yet");
     732                break;
     733        default:
     734                _error_("no case higher than 3 for SedimentlimitFlag");
     735        }
     736        /*Assign output pointer*/
     737        *ph_max=h_max;
     738}
     739/*}}}*/

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

    r17361 r17372  
    88/*Headers*/
    99#include "./Analysis.h"
    10 
     10class Node;
    1111class HydrologyDCEfficientAnalysis: public Analysis{
    1212
    … …  
    3434                IssmDouble EplSpecificStoring(Element* element);
    3535                IssmDouble SedimentStoring(Element* element);
    36                 IssmDouble GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss);
    3736                IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss);
    3837                IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss);
     38                void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask,Element* element);
     39                void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element);
     40                void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode);
    3941                void ComputeEPLThickness(FemModel* femmodel);
    4042};

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

    r17355 r17372  
    22#include "../toolkits/toolkits.h"
    33#include "../classes/classes.h"
     4#include "../classes/Node.h"
    45#include "../shared/shared.h"
    56#include "../modules/modules.h"
    … …  
    419420
    420421                for(int i=0;i<numnodes;i++){
    421                         basalelement->GetHydrologyDCInefficientHmax(&h_max,basalelement->GetNode(i));
     422                        GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->GetNode(i));
    422423                        if(values[i]>h_max) residual[i] = kappa*(values[i]-h_max);
    423424                        else                residual[i] = 0.;
    … …  
    496497        return h_max;
    497498}/*}}}*/
    498 
     499void  HydrologyDCInefficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/
     500       
     501        int        hmax_flag;
     502        IssmDouble h_max;
     503        IssmDouble rho_ice,rho_water;
     504        IssmDouble thickness,bed;
     505        /*Get the flag to the limitation method*/
     506        element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
     507       
     508        /*Switch between the different cases*/
     509        switch(hmax_flag){
     510        case 0:
     511                h_max=1.0e+10;
     512                break;
     513        case 1:
     514                element->FindParam(&h_max,HydrologydcSedimentlimitEnum);
     515                break;
     516        case 2:
     517                /*Compute max*/
     518                rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
     519                rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
     520                element->GetInputValue(&thickness,innode,ThicknessEnum);
     521                element->GetInputValue(&bed,innode,BedEnum);
     522                h_max=((rho_ice*thickness)/rho_water)+bed;
     523                break;
     524        case 3:
     525                _error_("Using normal stress  not supported yet");
     526                break;
     527        default:
     528                _error_("no case higher than 3 for SedimentlimitFlag");
     529        }
     530        /*Assign output pointer*/
     531        *ph_max=h_max;
     532}
     533/*}}}*/
    499534IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss){/*{{{*/
    500535

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

    r17349 r17372  
    88/*Headers*/
    99#include "./Analysis.h"
    10 
     10class Node;
    1111class HydrologyDCInefficientAnalysis: public Analysis{
    1212
    … …  
    3535                IssmDouble EplSpecificStoring(Element* element);
    3636                IssmDouble GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss);
     37                void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode);
    3738                IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss);
    3839                IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss);

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

    r17362 r17372  
    283283                virtual void UpdateConstraintsExtrudeFromTop(void)=0;
    284284
    285                 virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0;
    286                 virtual void HydrologyEPLGetMask(Vector<IssmDouble>* mask)=0;
    287                 virtual void HydrologyEPLGetActive(Vector<IssmDouble>* active)=0;
    288                 //              virtual void ComputeEPLThickness(void)=0;
    289 
    290285                virtual void   MigrateGroundingLine(IssmDouble* sheet_ungrounding)=0;
    291286                virtual void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding)=0;

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

    r17361 r17372  
    41554155/*}}}*/
    41564156#endif
    4157 
    4158 /*FUNCTION Penta::GetHydrologyDCInefficientHmax{{{*/
    4159 void  Penta::GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){
    4160 
    4161         if (!IsOnBed()) return;
    4162 
    4163         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    4164         tria->GetHydrologyDCInefficientHmax(ph_max,innode);
    4165         delete tria->material; delete tria;
    4166 }
    4167 /*}}}*/
    41684157/*FUNCTION Penta::GetSolutionFromInputsOneDof {{{*/
    41694158void Penta::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){
    … …  
    41974186}
    41984187/*}}}*/
    4199 /*FUNCTION Penta::HydrologyEPLGetActive {{{*/
    4200 void Penta::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){
    4201 
    4202         if (!IsOnBed()){
    4203                 return;
    4204         }
    4205         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    4206         tria->HydrologyEPLGetActive(active_vec);
    4207         delete tria->material; delete tria;
    4208 
    4209 }
    4210 /*}}}*/
    4211 /*FUNCTION Penta::HydrologyEPLGetMask{{{*/
    4212 void  Penta::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){
    4213 
    4214         if (!IsOnBed())return;
    4215 
    4216         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    4217         tria->HydrologyEPLGetMask(vec_mask);
    4218         delete tria->material; delete tria;
    4219 
    4220 }
    4221 /*}}}*/
    4222 /* /\*FUNCTION Penta::ComputeEPLThickness{{{*\/ */
    4223 /* void  Penta::ComputeEPLThickness(void){ */
    4224 
    4225 /*      int         i; */
    4226 /*      const int   numdof   = NDOF1 *NUMVERTICES; */
    4227 /*      const int   numdof2d = NDOF1 *NUMVERTICES2D; */
    4228 /*      bool        isefficientlayer; */
    4229 /*      IssmDouble  n,A,dt,init_thick; */
    4230 /*      IssmDouble  rho_water,rho_ice; */
    4231 /*      IssmDouble  gravity,latentheat,EPLgrad; */
    4232 /*      IssmDouble  EPL_N,epl_conductivity; */
    4233 /*      IssmDouble  activeEpl[numdof],thickness[numdof]; */
    4234 /*      IssmDouble  eplhead[numdof], old_thickness[numdof]; */
    4235 /*      IssmDouble  epl_slopeX[numdof],epl_slopeY[numdof]; */
    4236 /*      IssmDouble  ice_thickness[numdof],bed[numdof]; */
    4237 /*      Penta       *penta = NULL; */
    4238 /*      /\*If not on bed, return*\/ */
    4239 /*      if (!IsOnBed())return; */
    4240 
    4241 /*      /\*Get the flag to know if the efficient layer is present*\/ */
    4242 /*      this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); */
    4243 /*      this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); */
    4244 
    4245 /*      if(isefficientlayer){ */
    4246 /*              /\*For now, assuming just one way to compute EPL thickness*\/ */
    4247 /*              rho_water        = matpar->GetRhoWater(); */
    4248 /*              rho_ice          = matpar->GetRhoIce(); */
    4249 /*              gravity          = matpar->GetG(); */
    4250 /*              latentheat       = matpar->GetLatentHeat(); */
    4251 /*              epl_conductivity = matpar->GetEplConductivity(); */
    4252 /*              init_thick       = matpar->GetEplInitialThickness(); */
    4253 /*              n                = material->GetN(); */
    4254 /*              A                = material->GetA(); */
    4255                
    4256 /*              GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveNodeEnum); */
    4257 /*              GetInputListOnVertices(&eplhead[0],EplHeadEnum); */
    4258 /*              GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum);  */
    4259 /*              GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); */
    4260 /*              GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); */
    4261 /*              GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); */
    4262 /*              GetInputListOnVertices(&bed[0],BedEnum); */
    4263                
    4264 /*              for(int i=0;i<numdof2d;i++){ */
    4265 /*                      /\*Keeping thickness to 1 if EPL is not active*\/ */
    4266 /*                      if(activeEpl[i]==0.0){ */
    4267 /*                              thickness[i]=init_thick; */
    4268 /*                              thickness[i+numdof2d]=thickness[i]; */
    4269 /*                      } */
    4270 /*                      else{ */
    4271 
    4272 /*                              /\*Compute first the effective pressure in the EPL*\/ */
    4273 /*                              EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i]))); */
    4274 /*                              if(EPL_N<0.0)EPL_N=0.0; */
    4275 /*                              /\*Get then the gradient of EPL heads*\/ */
    4276 /*                              EPLgrad = epl_slopeX[i]+epl_slopeY[i]; */
    4277                                
    4278 /*                              /\*And proceed to the real thing*\/ */
    4279 /*                              thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))); */
    4280 /*                              thickness[i+numdof2d]=thickness[i]; */
    4281 /*                      } */
    4282 /*              } */
    4283 /*              penta=this; */
    4284 /*              for(;;){ */
    4285 
    4286 /*                      /\*Add input to the element: *\/                         */
    4287 /*                      penta->inputs->AddInput(new PentaInput(HydrologydcEplThicknessEnum,thickness,P1Enum)); */
    4288                        
    4289 /*                      /\*Stop if we have reached the surface*\/ */
    4290 /*                      if (penta->IsOnSurface()) break; */
    4291                        
    4292 /*                      /\* get upper Penta*\/ */
    4293 /*                      penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id); */
    4294 /*              } */
    4295 /*      } */
    4296 /* } */
    4297 /* /\*}}}*\/ */
    4298 
    42994188/*FUNCTION Penta::MigrateGroundingLine{{{*/
    43004189void  Penta::MigrateGroundingLine(IssmDouble* phi_ungrounding){

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

    r17361 r17372  
    241241                IssmDouble     StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa);
    242242
    243                 void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    244                 void           HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
    245                 void           HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
    246                 //              void           ComputeEPLThickness(void);
    247 
    248243                void           UpdateConstraintsExtrudeFromBase(void);
    249244                void           UpdateConstraintsExtrudeFromTop(void);

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

    r17361 r17372  
    150150                void        GetNormalFromLSF(IssmDouble *pnormal){_error_("not implemented yet");};
    151151
    152                 void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){_error_("not implemented yet");};
    153                 void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){_error_("not implemented yet");};
    154                 void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){_error_("not implemented yet");};
     152                //void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){_error_("not implemented yet");};
     153                //void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){_error_("not implemented yet");};
     154                //void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){_error_("not implemented yet");};
    155155                //              void    ComputeEPLThickness(void){_error_("not implemented yet");};
    156156               

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

    r17361 r17372  
    44254425}
    44264426/*}}}*/
    4427 /*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/
    4428 void  Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){
    4429        
    4430         int        hmax_flag;
    4431         IssmDouble h_max;
    4432         IssmDouble rho_ice,rho_water;
    4433         IssmDouble thickness,bed;
    4434         /*Get the flag to the limitation method*/
    4435         this->parameters->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
    4436        
    4437         /*Switch between the different cases*/
    4438         switch(hmax_flag){
    4439         case 0:
    4440                 h_max=1.0e+10;
    4441                 break;
    4442         case 1:
    4443                 parameters->FindParam(&h_max,HydrologydcSedimentlimitEnum);
    4444                 break;
    4445         case 2:
    4446                 rho_ice=matpar->GetRhoIce();
    4447                 rho_water=matpar->GetRhoFreshwater();
    4448                 this->GetInputValue(&thickness,innode,ThicknessEnum);
    4449                 this->GetInputValue(&bed,innode,BedEnum);
    4450                 h_max=((rho_ice*thickness)/rho_water)+bed;
    4451                 break;
    4452         case 3:
    4453                 _error_("Using normal stress  not supported yet");
    4454                 break;
    4455         default:
    4456                 _error_("no case higher than 3 for SedimentlimitFlag");
    4457         }
    4458         /*Assign output pointer*/
    4459         *ph_max=h_max;
    4460 }
    4461 /*}}}*/
    4462 /*FUNCTION Tria::HydrologyEPLGetActive {{{*/
    4463 void Tria::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){
    4464 
    4465         /*Constants*/
    4466         const int  numnodes = NUMVERTICES;
    4467         IssmDouble flag     = 0.;
    4468         IssmDouble active[numnodes];
    4469 
    4470         GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum);
    4471        
    4472         for(int i=0;i<numnodes;i++) flag+=active[i];
    4473 
    4474         if(flag>0.){
    4475                 for(int i=0;i<numnodes;i++){
    4476                         active_vec->SetValue(nodes[i]->Sid(),1.,INS_VAL);
    4477                 }
    4478         }
    4479         else{
    4480                 /*Do not do anything: at least one node is active for this element but this element is not solved for*/
    4481         }
    4482 }
    4483 /*}}}*/
    4484 /*FUNCTION Tria::HydrologyEPLGetMask{{{*/
    4485 void  Tria::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){
    4486 
    4487         /*Intermediaries*/
    4488         int         i,j;
    4489         const int   numdof         = NDOF1 *NUMVERTICES;
    4490         IssmDouble  init_thick;
    4491         IssmDouble  h_max;
    4492         IssmDouble  sedheadmin;
    4493         IssmDouble  epl_thickness[numdof];
    4494         IssmDouble  old_active[numdof];
    4495         IssmDouble  sedhead[numdof];
    4496         IssmDouble  eplhead[numdof];
    4497         IssmDouble  residual[numdof];
    4498 
    4499         bool       active_element;
    4500         Input* active_element_input=NULL;
    4501 
    4502         init_thick = matpar->GetEplInitialThickness();
    4503 
    4504         active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
    4505         active_element_input->GetInputValue(&active_element);
    4506 
    4507         GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum);       
    4508         GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum); 
    4509         GetInputListOnVertices(&sedhead[0],SedimentHeadEnum);
    4510         GetInputListOnVertices(&eplhead[0],EplHeadEnum);
    4511         GetInputListOnVertices(&residual[0],SedimentHeadResidualEnum);
    4512 
    4513         /*Get minimum sediment head of the element*/
    4514         sedheadmin=sedhead[0];
    4515         for(i=1;i<numdof;i++) if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i];
    4516 
    4517         for(i=0;i<numdof;i++){
    4518                 /*Activate EPL if residual is >0 */
    4519                 if(residual[i]>0.){
    4520                         vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL);
    4521                 }
    4522 
    4523                 /*If mask was already one, keep one*/
    4524                 else if(old_active[i]>0.){
    4525                         vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL);
    4526                         /*If epl thickness gets under , close the layer*/
    4527                         if(epl_thickness[i]<0.001*init_thick){
    4528                                 vec_mask->SetValue(nodes[i]->Sid(),0.,INS_VAL);
    4529                                 epl_thickness[i]=init_thick;
    4530                         }
    4531                 }
    4532                 /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
    4533                 this->GetHydrologyDCInefficientHmax(&h_max,this->nodes[i]);
    4534                 if(eplhead[i]>=h_max && active_element){
    4535                         for(j=0;j<numdof;j++){
    4536                                 if(old_active[j]>0.){
    4537                                         vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL);
    4538                                 }
    4539                                 /*Increase of the domain is on the downstream node in term of sediment head*/
    4540                                 if(sedhead[j] == sedheadmin){
    4541                                         vec_mask->SetValue(nodes[j]->Sid(),1.,INS_VAL);
    4542                                 }
    4543                         }
    4544                 }
    4545         }
    4546 }
    4547 /*}}}*/
    4548 //* *FUNCTION Tria::ComputeEPLThickness{{{*\/ */
    4549 /* void  Tria::ComputeEPLThickness(void){ */
    4550 
    4551 /*      const int   numdof         = NDOF1 *NUMVERTICES; */
    4552 /*      bool        isefficientlayer; */
    4553 /*      bool        active_element; */
    4554 /*      IssmDouble  n,A,dt,init_thick; */
    4555 /*      IssmDouble  rho_water,rho_ice; */
    4556 /*      IssmDouble  gravity,latentheat,EPLgrad2; */
    4557 /*      IssmDouble  EPL_N,epl_conductivity; */
    4558 /*      IssmDouble  thickness[numdof]; */
    4559 /*      IssmDouble  eplhead[numdof]; */
    4560 /*      IssmDouble  epl_slopeX[numdof],epl_slopeY[numdof]; */
    4561 /*      IssmDouble  old_thickness[numdof]; */
    4562 /*      IssmDouble  ice_thickness[numdof],bed[numdof]; */
    4563 
    4564 /*      Input* active_element_input=NULL; */
    4565 
    4566 /*      /\*Get the flag to know if the efficient layer is present*\/ */
    4567 /*      this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); */
    4568 /*      this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); */
    4569 
    4570 /*      if(isefficientlayer){ */
    4571 
    4572 /*              /\*For now, assuming just one way to compute EPL thickness*\/ */
    4573 /*              rho_water        = matpar->GetRhoWater(); */
    4574 /*              rho_ice          = matpar->GetRhoIce(); */
    4575 /*              gravity          = matpar->GetG(); */
    4576 /*              latentheat       = matpar->GetLatentHeat(); */
    4577 /*              epl_conductivity = matpar->GetEplConductivity(); */
    4578 /*              init_thick       = matpar->GetEplInitialThickness(); */
    4579 /*              n                = material->GetN(); */
    4580 /*              A                = material->GetAbar(); */
    4581                
    4582 /*              active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */
    4583 /*              active_element_input->GetInputValue(&active_element); */
    4584                        
    4585 /*              GetInputListOnVertices(&eplhead[0],EplHeadEnum); */
    4586 /*              GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum);  */
    4587 /*              GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); */
    4588 /*              GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); */
    4589 /*              GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); */
    4590 /*              GetInputListOnVertices(&bed[0],BedEnum); */
    4591                
    4592 /*              if(!active_element){ */
    4593 /*                      /\*Keeping thickness to initial value if EPL is not active*\/ */
    4594 /*                      for(int i=0;i<numdof;i++){ */
    4595 /*                              thickness[i]=init_thick; */
    4596 /*                      } */
    4597 /*              } */
    4598 /*              else{ */
    4599 /*                      for(int i=0;i<numdof;i++){ */
    4600                                
    4601 /*                              /\*Compute first the effective pressure in the EPL*\/ */
    4602 /*                              EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i]))); */
    4603 /*                              if(EPL_N<0.0)EPL_N=0.0; */
    4604 /*                              /\*Get then the square of the gradient of EPL heads*\/ */
    4605 /*                              EPLgrad2 = (epl_slopeX[i]*epl_slopeX[i]+epl_slopeY[i]*epl_slopeY[i]); */
    4606                                
    4607 /*                              /\*And proceed to the real thing*\/ */
    4608 /*                              thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))); */
    4609                                        
    4610 /*                              /\*Take care of otherthikening*\/ */
    4611 /*                              if(thickness[i]>10.0*init_thick){ */
    4612 /*                                      thickness[i] = 10.0*init_thick; */
    4613 /*                              } */
    4614 /*                      } */
    4615 /*              } */
    4616 /*              this->inputs->AddInput(new TriaInput(HydrologydcEplThicknessEnum,thickness,P1Enum)); */
    4617 /*      } */
    4618 /* } */
    4619 /*}}}*/
    46204427
    46214428#ifdef _HAVE_DAKOTA_

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

    r17361 r17372  
    253253
    254254                void           CreateHydrologyWaterVelocityInput(void);
    255                
    256                 void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    257                 void           HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
    258                 void           HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
    259                 //              void           ComputeEPLThickness(void);
    260255                /*}}}*/
    261256

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

    r17361 r17372  
    1515#include "./modules/modules.h"
    1616#include "../shared/Enum/Enum.h"
     17
     18#include "../analyses/analyses.h"
    1719
    1820/*module includes: {{{*/
    … …  
    13521354        Vector<IssmDouble>* active        = NULL;
    13531355        IssmDouble*         serial_active = NULL;
    1354 
    1355         /*Step 1: update maks, the mask might be extended by residual and/or using downstream sediment head*/
     1356       
     1357        HydrologyDCEfficientAnalysis* effanalysis =  new HydrologyDCEfficientAnalysis();
     1358        /*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/
    13561359        mask=new Vector<IssmDouble>(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
     1360
    13571361        for (int i=0;i<elements->Size();i++){
    13581362                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
    1359                 element->HydrologyEPLGetMask(mask);
     1363                effanalysis->HydrologyEPLGetMask(mask,element);
    13601364        }
    13611365
    … …  
    13731377        for (int i=0;i<elements->Size();i++){
    13741378                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
    1375                 element->HydrologyEPLGetActive(active);
     1379                effanalysis->HydrologyEPLGetActive(active,element);
    13761380        }
    13771381
    … …  
    13961400        }
    13971401        xDelete<IssmDouble>(serial_active);
     1402        delete effanalysis;
    13981403        int sum_counter;
    13991404        ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    … …  
    14071412}
    14081413/*}}}*/
    1409 /* void FemModel::HydrologyEPLThicknessx(void){ /\*{{{*\/ */
    1410 
    1411 /*      for (int i=0;i<elements->Size();i++){ */
    1412 /*              Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i)); */
    1413 /*              element->ComputeEPLThickness(); */
    1414 /*      } */
    1415 /* } */
    1416 /*}}}*/
    14171414void FemModel::UpdateConstraintsL2ProjectionEPLx(void){ /*{{{*/
    14181415
    14191416        Vector<IssmDouble>* active        = NULL;
    14201417        IssmDouble*         serial_active = NULL;
     1418        HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis();
    14211419
    14221420        /*update node activity. If one element is connected to mask=1, all nodes are active*/
    … …  
    14241422        for (int i=0;i<elements->Size();i++){
    14251423                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
    1426                 element->HydrologyEPLGetActive(active);
     1424                effanalysis->HydrologyEPLGetActive(active,element);
    14271425        }
    14281426
    … …  
    14311429        serial_active=active->ToMPISerial();
    14321430        delete active;
     1431        delete effanalysis;
    14331432
    14341433        /*Update node activation accordingly*/

  • issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp

    r17361 r17372  
    1313#include "../classes.h"
    1414#include "shared/shared.h"
     15#include "../../analyses/analyses.h"
    1516/*}}}*/
    1617
    … …  
    601602        IssmDouble h;
    602603        IssmDouble h_max;       
     604        HydrologyDCInefficientAnalysis* inefanalysis = NULL;
    603605
    604606        /*check that pengrid is not a clone (penalty to be added only once)*/
    … …  
    616618
    617619        /*Get sediment water head h*/
     620        inefanalysis = new HydrologyDCInefficientAnalysis();
    618621        element->GetInputValue(&h,node,SedimentHeadEnum);
    619         element->GetHydrologyDCInefficientHmax(&h_max,node);
     622        inefanalysis->GetHydrologyDCInefficientHmax(&h_max,element,node);
    620623        parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum);
    621624
    … …  
    646649
    647650        /*Assign output pointers:*/
     651        delete inefanalysis;
    648652        *punstable=unstable;
    649653}
    … …  
    651655/*FUNCTION Pengrid::PenaltyCreateKMatrixHydrologyDCInefficient {{{*/
    652656ElementMatrix* Pengrid::PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax){
    653 
    654657        IssmDouble    penalty_factor;
    655658
    … …  
    672675        IssmDouble h_max;
    673676        IssmDouble penalty_factor;
     677        HydrologyDCInefficientAnalysis* inefanalysis = NULL;
    674678
    675679        /*Initialize Element matrix and return if necessary*/
    676680        if(!this->active) return NULL;
    677681        ElementVector* pe=new ElementVector(&node,1,this->parameters);
     682        inefanalysis = new HydrologyDCInefficientAnalysis();
    678683
    679684        /*Retrieve parameters*/
    … …  
    681686
    682687        /*Get h_max and compute penalty*/
    683         element->GetHydrologyDCInefficientHmax(&h_max,node);
     688        inefanalysis->GetHydrologyDCInefficientHmax(&h_max,element,node);
     689
    684690        pe->values[0]=kmax*pow(10.,penalty_factor)*h_max;
    685691
    686692        /*Clean up and return*/
     693        delete inefanalysis;
    687694        return pe;
    688695}

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

    r17362 r17372  
    156156                }
    157157                /* }}} *//*End of the global sediment loop*/
    158 
    159158                /* {{{ *//*Now dealing with the EPL in the same way*/
    160159                if(isefficientlayer){
    … …  
    170169                                ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter);
    171170                                ug_epl->Copy(ug_epl_sub_iter);
    172 
    173171                                /* {{{ *//*Retrieve the EPL head slopes and compute EPL Thickness*/
    174172                                if(VerboseSolution()) _printf0_("computing EPL Head slope...\n");

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

    r16716 r17372  
    2424        femmodel->UpdateConstraintsx();
    2525        SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
    26 
    2726        CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
    2827        Reduceloadx(pf, Kfs, ys); delete Kfs;
Note: See TracChangeset for help on using the changeset viewer.