Changeset 16685


Ignore:
Timestamp:
11/08/13 11:50:50 (11 years ago)
Author:
bdef
Message:

Adding stuff to compute Epl slopes only on active EPL

Location:
issm/trunk-jpl/src/c
Files:
2 added
13 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r16675 r16685  
    465465                                                        ./analyses/HydrologyShreveAnalysis.h\
    466466                                                        ./analyses/HydrologyShreveAnalysis.cpp\
     467                                                        ./analyses/L2ProjectionEPLAnalysis.h\
     468                                                ./analyses/L2ProjectionEPLAnalysis.cpp\
    467469                                                        ./cores/hydrology_core.cpp\
    468470                                                        ./solutionsequences/solutionsequence_hydro_nonlinear.cpp

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

    r16534 r16685  
    2020                case BalancevelocityAnalysisEnum : return new BalancevelocityAnalysis();
    2121                case L2ProjectionBaseAnalysisEnum : return new L2ProjectionBaseAnalysis();
     22                case L2ProjectionEPLAnalysisEnum : return new L2ProjectionEPLAnalysis();
    2223                case DamageEvolutionAnalysisEnum : return new DamageEvolutionAnalysis();
    2324                case StressbalanceAnalysisEnum : return new StressbalanceAnalysis();

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

    r16534 r16685  
    3232#include "./StressbalanceVerticalAnalysis.h"
    3333#include "./L2ProjectionBaseAnalysis.h"
     34#include "./L2ProjectionEPLAnalysis.h"
    3435#include "./ThermalAnalysis.h"
    35 
    36 #include "EnumToAnalysis.h"
     36#include "./EnumToAnalysis.h"
    3737#endif

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

    r16681 r16685  
    154154                virtual void HydrologyEPLGetActive(Vector<IssmDouble>* active)=0;
    155155                virtual void ComputeEPLThickness(void)=0;
     156                virtual void UpdateConstraintsL2ProjectionEPL(void)=0;
    156157                #endif
    157158

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

    r16681 r16685  
    481481          case HydrologyDCEfficientAnalysisEnum:
    482482                        return CreateKMatrixHydrologyDCEfficient();
     483                        break;
     484          case L2ProjectionEPLAnalysisEnum:
     485                        return CreateEPLDomainMassMatrix();
    483486                        break;
    484487                #endif
    … …  
    697700                        return CreatePVectorHydrologyDCEfficient();
    698701                        break;
     702          case L2ProjectionEPLAnalysisEnum:
     703                        return CreatePVectorL2ProjectionEPL();
     704                        break;
    699705                #endif
    700706                default:
    … …  
    22582264        case HydrologyDCEfficientAnalysisEnum:
    22592265                InputUpdateFromSolutionOneDofCollapsed(solution,EplHeadEnum);
     2266                break;
     2267        case L2ProjectionEPLAnalysisEnum:
     2268                this->parameters->FindParam(&inputenum,InputToL2ProjectEnum);
     2269                InputUpdateFromSolutionOneDofCollapsed(solution,inputenum);
    22602270                break;
    22612271        #endif
    … …  
    1062210632}
    1062310633/*}}}*/
     10634/*FUNCTION Penta::CreateEPLDomainMassMatrix {{{*/
     10635ElementMatrix* Penta::CreateEPLDomainMassMatrix(void){
     10636
     10637        if (!IsOnBed()) return NULL;
     10638
     10639        Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
     10640        ElementMatrix* Ke=tria->CreateEPLDomainMassMatrix();
     10641        delete tria->material; delete tria;
     10642
     10643        /*clean up and return*/
     10644        return Ke;
     10645}
     10646/*}}}*/
    1062410647/*FUNCTION Penta::CreatePVectorHydrologyDCInefficient {{{*/
    1062510648ElementVector* Penta::CreatePVectorHydrologyDCInefficient(void){
    … …  
    1064410667        Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    1064510668        ElementVector* pe=tria->CreatePVectorHydrologyDCEfficient();
     10669        delete tria->material; delete tria;
     10670
     10671        /*Clean up and return*/
     10672        return pe;
     10673}
     10674/*}}}*/
     10675/*FUNCTION Penta::CreatePVectorL@ProjectionEPL {{{*/
     10676ElementVector* Penta::CreatePVectorL2ProjectionEPL(void){
     10677
     10678        if (!IsOnBed()) return NULL;
     10679
     10680        /*Call Tria function*/
     10681        Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
     10682        ElementVector* pe=tria->CreatePVectorL2ProjectionEPL();
    1064610683        delete tria->material; delete tria;
    1064710684
    … …  
    1077510812                                /*Compute first the effective pressure in the EPL*/
    1077610813                                EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
    10777                                
     10814                                if(EPL_N<0.0)EPL_N=0.0;
    1077810815                                /*Get then the gradient of EPL heads*/
    1077910816                                EPLgrad = epl_slopeX[i]+epl_slopeY[i];
    1078010817                               
    1078110818                                /*And proceed to the real thing*/
    10782                                 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)))*EPL_N);
     10819                                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)));
    1078310820                                thickness[i+numdof2d]=thickness[i];
     10821                                printf("N, %e - thick term2, %e \n",EPL_N,-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)));
     10822                                printf("old thick, %e - thick , %e \n",old_thickness[i],thickness[i]);
    1078410823                        }
    1078510824                }
    … …  
    1086610905        /*Free ressources:*/
    1086710906        xDelete<int>(doflist);
     10907}
     10908/*}}}*/
     10909/*FUNCTION Penta::UpdateConstraintsL2ProjectionEPL{{{*/
     10910void  Penta::UpdateConstraintsL2ProjectionEPL(void){
     10911
     10912        IssmDouble activeEpl[NUMVERTICES];
     10913
     10914
     10915        if(!IsOnBed()){
     10916                for(int i=0;i<this->NumberofNodes();i++)this->nodes[i]->Deactivate();
     10917        }
     10918        else{
     10919                GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
     10920                for(int i=0;i<3;i++){
     10921                        if(!activeEpl[i])this->nodes[i]->Deactivate();
     10922                }
     10923        }
    1086810924}
    1086910925/*}}}*/

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

    r16681 r16685  
    325325                ElementMatrix* CreateKMatrixHydrologyDCInefficient(void);
    326326                ElementMatrix* CreateKMatrixHydrologyDCEfficient(void);
     327                ElementMatrix* CreateEPLDomainMassMatrix(void);
    327328                ElementVector* CreatePVectorHydrologyDCInefficient(void);
    328329                ElementVector* CreatePVectorHydrologyDCEfficient(void);
    329                 void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    330                 void    GetHydrologyTransfer(Vector<IssmDouble>* transfer);
    331                 void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
    332                 void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
    333                 void    InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
    334                 void    ComputeEPLThickness(void);
     330                ElementVector* CreatePVectorL2ProjectionEPL(void);
     331                void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
     332                void           GetHydrologyTransfer(Vector<IssmDouble>* transfer);
     333                void           HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
     334                void           HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
     335                void           InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
     336                void           ComputeEPLThickness(void);
     337                void           UpdateConstraintsL2ProjectionEPL(void);
    335338                #endif
    336339

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

    r16681 r16685  
    117117                void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){_error_("not implemented yet");};
    118118                void    ComputeEPLThickness(void){_error_("not implemented yet");};
     119                void    UpdateConstraintsL2ProjectionEPL(){_error_("not implemented");};
    119120                #endif
    120121                void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type){_error_("not implemented yet");};

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

    r16681 r16685  
    288288                case HydrologyDCEfficientAnalysisEnum:
    289289                        return CreateKMatrixHydrologyDCEfficient();
     290                        break;
     291          case L2ProjectionEPLAnalysisEnum:
     292                        return CreateEPLDomainMassMatrix();
    290293                        break;
    291294                #endif
    … …  
    513516                        return CreatePVectorHydrologyDCEfficient();
    514517                        break;
     518          case L2ProjectionEPLAnalysisEnum:
     519                        return CreatePVectorL2ProjectionEPL();
     520                        break;
    515521                #endif
    516522                #ifdef _HAVE_BALANCED_
    … …  
    586592                case BedSlopeXEnum:     input2 = inputs->GetInput(BedEnum);     _assert_(input2); break;
    587593                case BedSlopeYEnum:     input2 = inputs->GetInput(BedEnum);     _assert_(input2); break;
    588                 case EplHeadSlopeXEnum: input2 = inputs->GetInput(EplHeadEnum); _assert_(input2); break;
    589                 case EplHeadSlopeYEnum: input2 = inputs->GetInput(EplHeadEnum); _assert_(input2); break;
    590594                default: input = inputs->GetInput(input_enum);
    591595        }
    … …  
    602606                if(input2) input2->GetInputDerivativeValue(&slopes[0],&xyz_list[0][0],gauss);
    603607                switch(input_enum){
    604                         case SurfaceSlopeXEnum: case BedSlopeXEnum: case EplHeadSlopeXEnum: value = slopes[0]; break;
    605                         case SurfaceSlopeYEnum: case BedSlopeYEnum: case EplHeadSlopeYEnum: value = slopes[1]; break;
     608                        case SurfaceSlopeXEnum: case BedSlopeXEnum: value = slopes[0]; break;
     609                        case SurfaceSlopeYEnum: case BedSlopeYEnum: value = slopes[1]; break;
    606610                        default: input->GetInputValue(&value,gauss);
    607611                }
    … …  
    17841788                        InputUpdateFromSolutionOneDof(solution,EplHeadEnum);
    17851789                        break;
     1790                case L2ProjectionEPLAnalysisEnum:
     1791                        this->parameters->FindParam(&extrusioninput,InputToL2ProjectEnum);
     1792                        InputUpdateFromSolutionOneDof(solution,extrusioninput);
     1793                        break;
    17861794                #endif
    17871795                #ifdef _HAVE_DAMAGE_
    … …  
    29012909}
    29022910/*}}}*/
    2903 
     2911/*FUNCTION Tria::UpdateConstraintsL2ProjectionEPL{{{*/
     2912void  Tria::UpdateConstraintsL2ProjectionEPL(void){
     2913
     2914        IssmDouble activeEpl[NUMVERTICES];
     2915 
     2916        GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
     2917        for(int i=0;i<NUMVERTICES;i++){
     2918                if(!activeEpl[i])this->nodes[i]->Deactivate();
     2919        }
     2920}
     2921/*}}}*/
    29042922#ifdef _HAVE_RESPONSES_
    29052923/*FUNCTION Tria::AverageOntoPartition {{{*/
    … …  
    66836701}
    66846702/*}}}*/
     6703
     6704/*FUNCTION Tria::CreatEPLDomainMassMatrix {{{*/
     6705ElementMatrix* Tria::CreateEPLDomainMassMatrix(void){
     6706
     6707        /* Intermediaries */
     6708        IssmDouble  D,Jdet;
     6709        IssmDouble  xyz_list[NUMVERTICES][3];
     6710
     6711        /*Fetch number of nodes and dof for this finite element*/
     6712        int numnodes = this->NumberofNodes();
     6713
     6714        /*Initialize Element matrix and vectors*/
     6715        ElementMatrix* Ke    = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     6716        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     6717
     6718        /*Retrieve all inputs and parameters*/
     6719        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6720
     6721        /* Start looping on the number of gaussian points: */
     6722        GaussTria* gauss=new GaussTria(2);
     6723        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6724
     6725                gauss->GaussPoint(ig);
     6726
     6727                GetNodalFunctions(basis,gauss);
     6728                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     6729                D=gauss->weight*Jdet;
     6730
     6731                TripleMultiply(basis,1,numnodes,1,
     6732                                        &D,1,1,0,
     6733                                        basis,1,numnodes,0,
     6734                                        &Ke->values[0],1);
     6735        }
     6736
     6737        /*Clean up and return*/
     6738        delete gauss;
     6739        xDelete<IssmDouble>(basis);
     6740        return Ke;
     6741}
     6742/*}}}*/
    66856743/*FUNCTION Tria::CreatePVectorHydrologyShreve {{{*/
    66866744ElementVector* Tria::CreatePVectorHydrologyShreve(void){
    … …  
    67086766        Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
    67096767        Input* old_watercolumn_input=inputs->GetInput(WaterColumnOldEnum);         _assert_(old_watercolumn_input);
    6710 
    67116768        /*Initialize basal_melting_correction_g to 0, do not forget!:*/
    67126769        /* Start  looping on the number of gaussian points: */
    … …  
    68646921                residual_input->GetInputValue(&residual,gauss);
    68656922                pe->values[iv]+=residual/connectivity;
     6923        }
     6924
     6925        /*Clean up and return*/
     6926        xDelete<IssmDouble>(basis);
     6927        delete gauss;
     6928        return pe;
     6929}
     6930/*}}}*/
     6931/*FUNCTION Tria::CreatePVectorL2ProjectionEPL {{{*/
     6932ElementVector* Tria::CreatePVectorL2ProjectionEPL(void){
     6933
     6934        /*Intermediaries */
     6935        int        i,input_enum;
     6936        IssmDouble Jdet,value;
     6937        IssmDouble xyz_list[NUMVERTICES][3];
     6938        IssmDouble slopes[2];
     6939        Input*     input  = NULL;
     6940        Input*     input2 = NULL;
     6941
     6942        /*Fetch number of nodes and dof for this finite element*/
     6943        int numnodes = this->NumberofNodes();
     6944
     6945        /*Initialize Element vector*/
     6946        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     6947        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     6948
     6949        /*Retrieve all inputs and parameters*/
     6950        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6951        this->parameters->FindParam(&input_enum,InputToL2ProjectEnum);
     6952        switch(input_enum){
     6953                case EplHeadSlopeXEnum: input2 = inputs->GetInput(SurfaceEnum); _assert_(input2); break;
     6954                case EplHeadSlopeYEnum: input2 = inputs->GetInput(SurfaceEnum); _assert_(input2); break;
     6955        default: input = inputs->GetInput(input_enum);
     6956        }
     6957
     6958        /* Start  looping on the number of gaussian points: */
     6959        GaussTria* gauss=new GaussTria(2);
     6960        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6961
     6962                gauss->GaussPoint(ig);
     6963
     6964                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     6965                GetNodalFunctions(basis,gauss);
     6966
     6967                if(input2) input2->GetInputDerivativeValue(&slopes[0],&xyz_list[0][0],gauss);
     6968                switch(input_enum){
     6969                        case EplHeadSlopeXEnum: value = slopes[0]; break;
     6970                        case EplHeadSlopeYEnum: value = slopes[1]; break;
     6971                        default: input->GetInputValue(&value,gauss);
     6972                }
     6973
     6974                for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*value*basis[i];
    68666975        }
    68676976
    … …  
    71277236        IssmDouble  h_max;
    71287237        IssmDouble  sedheadmin;
     7238        IssmDouble  epl_thickness[numdof];
    71297239        IssmDouble  old_active[numdof];
    71307240        IssmDouble  sedhead[numdof];
    … …  
    71337243
    71347244        GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveEnum);   
     7245        GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum); 
    71357246        GetInputListOnVertices(&sedhead[0],SedimentHeadEnum);
    71367247        GetInputListOnVertices(&eplhead[0],EplHeadEnum);
    … …  
    71497260                /*If mask was alread one, keep one*/
    71507261                else if(old_active[i]>0.){
    7151                         vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL);
    7152                 }
    7153 
     7262                        if(epl_thickness[i]>0.0){
     7263                                vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL);
     7264                        }
     7265                        else{
     7266                                vec_mask->SetValue(nodes[i]->Sid(),0.,INS_VAL);
     7267                        }
     7268                }
    71547269                /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
    71557270                this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
    … …  
    72177332                                /*Compute first the effective pressure in the EPL*/
    72187333                                EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
    7219                                
     7334                                if(EPL_N<0.0)EPL_N=0.0;
    72207335                                /*Get then the gradient of EPL heads*/
    72217336                                EPLgrad = epl_slopeX[i]+epl_slopeY[i];
    72227337                               
    72237338                                /*And proceed to the real thing*/
    7224                                 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)))*EPL_N);
    7225 
     7339                                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)));
     7340                                if(this->nodes[i]->id==1553)printf("term1  %e, term2 %e\n",+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*pow(EPLgrad,2.0), -2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)));
    72267341                        }
    72277342                }

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

    r16681 r16685  
    301301                ElementMatrix* CreateKMatrixHydrologyDCInefficient(void);
    302302                ElementMatrix* CreateKMatrixHydrologyDCEfficient(void);
     303                ElementMatrix* CreateEPLDomainMassMatrix(void);
    303304                ElementVector* CreatePVectorHydrologyShreve(void);
    304305                ElementVector* CreatePVectorHydrologyDCInefficient(void);
    305306                ElementVector* CreatePVectorHydrologyDCEfficient(void);
    306                 void    CreateHydrologyWaterVelocityInput(void);
    307                 void      InputUpdateFromSolutionHydrology(IssmDouble* solution);
    308                 void      InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution);
    309                 void    InputUpdateFromSolutionHydrologyDC(IssmDouble* solution);
    310                 void      InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
    311                 void      InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution);
    312                 void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    313                 void    GetHydrologyTransfer(Vector<IssmDouble>* transfer);
    314                 void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
    315                 void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
    316                 void    ComputeEPLThickness(void);
    317                 bool    AllActive(void);
    318                 bool    AnyActive(void);
     307                ElementVector* CreatePVectorL2ProjectionEPL(void);
     308                void           CreateHydrologyWaterVelocityInput(void);
     309                void             InputUpdateFromSolutionHydrology(IssmDouble* solution);
     310                void             InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution);
     311                void           InputUpdateFromSolutionHydrologyDC(IssmDouble* solution);
     312                void             InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
     313                void             InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution);
     314                void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
     315                void           GetHydrologyTransfer(Vector<IssmDouble>* transfer);
     316                void           HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
     317                void           HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
     318                void           ComputeEPLThickness(void);
     319                void           UpdateConstraintsL2ProjectionEPL(void);
     320                bool           AllActive(void);
     321                bool           AnyActive(void);
    319322                #endif
    320323

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

    r16600 r16685  
    261261         * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the
    262262         * Slope configuration.*/
    263 
    264263        int found=-1;
    265264        for(int i=0;i<nummodels;i++){
     265       
    266266                if (analysis_type_list[i]==configuration_type){
    267267                        found=i;
    … …  
    14201420}
    14211421/*}}}*/
    1422 #endif
     1422
     1423void FemModel::UpdateConstraintsL2ProjectionEPLx(void){ /*{{{*/
     1424
     1425        for(int i=0;i<elements->Size();i++){
     1426                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     1427                element->UpdateConstraintsL2ProjectionEPL();
     1428        }
     1429
     1430}
     1431/*}}}*/#endif

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

    r16600 r16685  
    103103                void HydrologyEPLupdateDomainx(void);
    104104                void HydrologyEPLThicknessx(void);
     105                void UpdateConstraintsL2ProjectionEPLx(void);
    105106                #endif
    106107};

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

    r16612 r16685  
    5656
    5757                case HydrologySolutionEnum:
    58                         numanalyses=4;
     58                        numanalyses=5;
    5959                        analyses=xNew<int>(numanalyses);
    6060                        analyses[0]=HydrologyShreveAnalysisEnum;
    … …  
    6262                        analyses[2]=HydrologyDCEfficientAnalysisEnum;
    6363                        analyses[3]=L2ProjectionBaseAnalysisEnum;
     64                        analyses[4]=L2ProjectionEPLAnalysisEnum;
    6465                        break;
    6566

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

    r16638 r16685  
    4646        femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum);
    4747        femmodel->parameters->FindParam(&time,TimeEnum);
     48        /*FIXME, hardcoded, put on an enum*/
    4849        hydro_maxiter=150;
    4950        hydrocount=1;
    … …  
    9899                                if(num_unstable_constraints==0) sedconverged = true;
    99100                                if (sedcount>=hydro_maxiter){
     101                                        /*Hacking to get the results of non converged runs*/
     102                                        //                                      sedconverged = true;
    100103                                        _error_("   maximum number of Sediment iterations (" << hydro_maxiter << ") exceeded");
    101104                                }
    … …  
    116119
    117120                        /*Start by retrieving the EPL head slopes*/
    118                         if(VerboseSolution()) _printf0_("computing EPL Head slope...\n");
    119                         femmodel->SetCurrentConfiguration(L2ProjectionBaseAnalysisEnum);
    120                         femmodel->parameters->SetParam(EplHeadSlopeXEnum,InputToL2ProjectEnum);
    121                         solutionsequence_linear(femmodel);
    122                         femmodel->parameters->SetParam(EplHeadSlopeYEnum,InputToL2ProjectEnum);
    123                         solutionsequence_linear(femmodel);
     121                        /* if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); */
     122                        /* femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum); */
     123                        /* femmodel->parameters->SetParam(EplHeadSlopeXEnum,InputToL2ProjectEnum); */
     124                        /* solutionsequence_linear(femmodel); */
     125                        /* femmodel->parameters->SetParam(EplHeadSlopeYEnum,InputToL2ProjectEnum); */
     126                        /* solutionsequence_linear(femmodel); */
    124127                       
    125128                        femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
    … …  
    134137                        eplconverged = false;
    135138                        for(;;){
     139
     140                        /*Start by retrieving the EPL head slopes*/
     141                                if(VerboseSolution()) _printf0_("computing EPL Head slope...\n");
     142                                femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum);
     143                                femmodel->UpdateConstraintsL2ProjectionEPLx();
     144                                femmodel->parameters->SetParam(EplHeadSlopeXEnum,InputToL2ProjectEnum);
     145                                solutionsequence_linear(femmodel);
     146                                femmodel->parameters->SetParam(EplHeadSlopeYEnum,InputToL2ProjectEnum);
     147                                solutionsequence_linear(femmodel);
     148                                femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
     149
    136150                                femmodel->HydrologyEPLThicknessx();
     151                                //updating mask after the computation of the epl thickness
     152                                femmodel->HydrologyEPLupdateDomainx();
    137153                                femmodel->HydrologyTransferx();
    138154                                SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
    … …  
    155171                                        if(num_unstable_constraints==0) eplconverged = true;
    156172                                        if (eplcount>=hydro_maxiter){
     173                                        /*Hacking to get the results of non converged runs*/
     174                                        //eplconverged = true;
    157175                                                _error_("   maximum number of EPL iterations (" << hydro_maxiter << ") exceeded");
    158176                                        }
    … …  
    214232                        else _printf0_(setw(50) << left << "   Convergence criterion:" << ndu_sed/nu_sed*100 << " %\n");
    215233                        if (hydrocount>=hydro_maxiter){
    216                                 _error_("   maximum number for hydrological global iterations (" << hydro_maxiter << ") exceeded");
     234                                        /*Hacking to get the results of non converged runs*/
     235                                        //hydroconverged = true;
     236                                        _error_("   maximum number for hydrological global iterations (" << hydro_maxiter << ") exceeded");
    217237                        }
    218238                }
Note: See TracChangeset for help on using the changeset viewer.