Changeset 14605


Ignore:
Timestamp:
04/16/13 15:10:47 (12 years ago)
Author:
seroussi
Message:

NEW: added subelement migration implementation in 3d

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

Legend:

Unmodified
Added
Removed

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

    r14591 r14605  
    826826        /*Assign output pointers:*/
    827827        *pdoflist=doflist;
     828}
     829/*}}}*/
     830/*FUNCTION Penta::GetGroundedPortion{{{*/
     831IssmDouble Penta::GetGroundedPortion(IssmDouble* xyz_list){
     832        /*Computeportion of the element that is grounded*/
     833
     834        bool               mainlyfloating = true;
     835        const IssmPDouble  epsilon= 1.e-15;
     836        IssmDouble         phi,s1,s2,area_init,area_grounded;
     837        IssmDouble         gl[3];
     838        IssmDouble         xyz_bis[3][3];
     839
     840        /*Recover parameters and values*/
     841        GetInputListOnVertices(&gl[0],GLlevelsetEnum);
     842       
     843        /*Be sure that values are not zero*/
     844        if(gl[0]==0) gl[0]=gl[0]+epsilon;
     845        if(gl[1]==0) gl[1]=gl[1]+epsilon;
     846        if(gl[2]==0) gl[2]=gl[2]+epsilon;
     847
     848        /*Check that not all nodes are grounded or floating*/
     849        if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
     850                phi=1;
     851        }
     852        else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
     853                phi=0;
     854        }
     855        else{
     856                /*Figure out if two nodes are floating or grounded*/
     857                if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=false;
     858
     859                if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     860                        /*Coordinates of point 2: same as initial point 2*/
     861                        xyz_bis[2][0]=*(xyz_list+3*2+0);
     862                        xyz_bis[2][1]=*(xyz_list+3*2+1);
     863                        xyz_bis[2][2]=*(xyz_list+3*2+2);
     864
     865                        /*Portion of the segments*/
     866                        s1=gl[2]/(gl[2]-gl[1]);
     867                        s2=gl[2]/(gl[2]-gl[0]);
     868
     869                        /*New point 1*/
     870                        xyz_bis[1][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*2+0));
     871                        xyz_bis[1][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*2+1));
     872                        xyz_bis[1][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*2+2));
     873
     874                        /*New point 0*/
     875                        xyz_bis[0][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*2+0));
     876                        xyz_bis[0][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*2+1));
     877                        xyz_bis[0][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*2+2));
     878                }
     879                else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
     880                        /*Coordinates of point 0: same as initial point 2*/
     881                        xyz_bis[0][0]=*(xyz_list+3*0+0);
     882                        xyz_bis[0][1]=*(xyz_list+3*0+1);
     883                        xyz_bis[0][2]=*(xyz_list+3*0+2);
     884
     885                        /*Portion of the segments*/
     886                        s1=gl[0]/(gl[0]-gl[1]);
     887                        s2=gl[0]/(gl[0]-gl[2]);
     888
     889                        /*New point 1*/
     890                        xyz_bis[1][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*0+0));
     891                        xyz_bis[1][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*0+1));
     892                        xyz_bis[1][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*0+2));
     893
     894                        /*New point 2*/
     895                        xyz_bis[2][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*2+0)-*(xyz_list+3*0+0));
     896                        xyz_bis[2][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*2+1)-*(xyz_list+3*0+1));
     897                        xyz_bis[2][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*2+2)-*(xyz_list+3*0+2));
     898                }
     899                else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
     900                        /*Coordinates of point 1: same as initial point 2*/
     901                        xyz_bis[1][0]=*(xyz_list+3*1+0);
     902                        xyz_bis[1][1]=*(xyz_list+3*1+1);
     903                        xyz_bis[1][2]=*(xyz_list+3*1+2);
     904
     905                        /*Portion of the segments*/
     906                        s1=gl[1]/(gl[1]-gl[0]);
     907                        s2=gl[1]/(gl[1]-gl[2]);
     908
     909                        /*New point 0*/
     910                        xyz_bis[0][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*1+0));
     911                        xyz_bis[0][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*1+1));
     912                        xyz_bis[0][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*1+2));
     913
     914                        /*New point 2*/
     915                        xyz_bis[2][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*2+0)-*(xyz_list+3*1+0));
     916                        xyz_bis[2][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*2+1)-*(xyz_list+3*1+1));
     917                        xyz_bis[2][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*2+2)-*(xyz_list+3*1+2));
     918                }
     919
     920                /*Compute fraction of grounded element*/
     921                GetJacobianDeterminant(&area_init, xyz_list,NULL);
     922                GetJacobianDeterminant(&area_grounded, &xyz_bis[0][0],NULL);
     923                if(mainlyfloating==true) area_grounded=area_init-area_grounded;
     924                phi=area_grounded/area_init;
     925        }
     926
     927        if(phi>1 || phi<0) _error_("Error. Problem with portion of grounded element: value should be between 0 and 1");
     928
     929        return phi;
    828930}
    829931/*}}}*/
    … …  
    20532155                                name==WaterfractionEnum||
    20542156                                name==FrictionCoefficientEnum ||
     2157                                name==GLlevelsetEnum ||
    20552158                                name==GradientEnum ||
    20562159                                name==OldGradientEnum  ||
    … …  
    21812284        int     i,migration_style;
    21822285        bool    floatingelement = false;
     2286        bool    groundedelement = false;
    21832287        IssmDouble  bed_hydro,yts,gl_melting_rate;
    21842288        IssmDouble  rho_water,rho_ice,density;
    2185         IssmDouble  melting[NUMVERTICES];
     2289        IssmDouble  melting[NUMVERTICES],phi[NUMVERTICES];
    21862290        IssmDouble  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
    21872291
    … …  
    21962300        GetInputListOnVertices(&b[0],BedEnum);
    21972301        GetInputListOnVertices(&ba[0],BathymetryEnum);
     2302        if(migration_style==SubelementMigrationEnum) GetInputListOnVertices(&phi[0],GLlevelsetEnum);
    21982303        rho_water=matpar->GetRhoWater();
    21992304        rho_ice=matpar->GetRhoIce();
    … …  
    22172322                        if (bed_hydro>ba[i]){
    22182323                                /*Unground only if the element is connected to the ice shelf*/
    2219                                 if(migration_style==AgressiveMigrationEnum){
     2324                                if(migration_style==AgressiveMigrationEnum || migration_style==SubelementMigrationEnum){
    22202325                                        s[i]=(1-density)*h[i];
    22212326                                        b[i]=-density*h[i];
    … …  
    22292334                                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
    22302335                                }
     2336                                else{
     2337                                        if(migration_style!=SoftMigrationEnum) _error_("Error: migration should be Aggressive, Soft or Subelement");
     2338                                }
    22312339                        }
    22322340                }
    22332341        }
    22342342
    2235         /*If at least one vertex is now floating, the element is now floating*/
    2236         for(i=0;i<NUMVERTICES;i++){
    2237                 if(nodes[i]->IsFloating()){
    2238                         floatingelement=true;
    2239                         break;
     2343        /*SubelementMigrationEnum: if one grounded, all grounded*/
     2344        if(migration_style==SubelementMigrationEnum){
     2345                for(i=0;i<NUMVERTICES;i++){
     2346                        if(nodes[i]->IsGrounded()){
     2347                                groundedelement=true;
     2348                                break;
     2349                        }
     2350                }
     2351                floatingelement=!groundedelement;
     2352        }
     2353        else{
     2354                for(i=0;i<NUMVERTICES;i++){
     2355                        if(nodes[i]->IsFloating()){
     2356                                floatingelement=true;
     2357                                break;
     2358                        }
    22402359                }
    22412360        }
    … …  
    22512370        this->inputs->AddInput(new PentaP1Input(BedEnum,&b[0]));
    22522371   this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,floatingelement));
     2372
     2373        /*Recalculate phi*/
     2374        if(migration_style==SubelementMigrationEnum){
     2375                for(i=0;i<NUMVERTICES;i++) phi[i]=h[i]+ba[i]/density;
     2376                this->inputs->AddInput(new PentaP1Input(GLlevelsetEnum,&phi[0]));
     2377                this->InputExtrude(GLlevelsetEnum,ElementEnum);
     2378        }
    22532379
    22542380        /*Extrude inputs*/
    … …  
    67676893        /*Intermediaries */
    67686894        int       i,j;
    6769         int       analysis_type;
     6895        int       analysis_type,migration_style;
    67706896        IssmDouble xyz_list[NUMVERTICES][3];
    67716897        IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0};
    67726898        IssmDouble slope_magnitude,alpha2,Jdet;
     6899        IssmDouble phi=1.0;
    67736900        IssmDouble slope[3]={0.0,0.0,0.0};
    67746901        IssmDouble MAXSLOPE=.06; // 6 %
    … …  
    67976924        friction=new Friction("2d",inputs,matpar,analysis_type);
    67986925
     6926        /*Recover portion of element that is grounded*/
     6927        if(migration_style==SubelementMigrationEnum) phi=this->GetGroundedPortion(&xyz_list[0][0]);
     6928
    67996929        /* Start  looping on the number of gaussian points: */
    68006930        gauss=new GaussPenta(0,1,2,2);
    … …  
    68096939                friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
    68106940                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
     6941                if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
    68116942
    68126943                // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,

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

    r14589 r14605  
    181181                ElementVector* CreatePVectorPrognostic(void);
    182182                ElementVector* CreatePVectorSlope(void);
    183                 void      GetDofList(int** pdoflist,int approximation_enum,int setenum);
    184                 void      GetVertexPidList(int* doflist);
    185                 void    GetVertexSidList(int* sidlist);
    186                 void    GetConnectivityList(int* connectivity);
    187                 int     GetElementType(void);
    188                 void    GetElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
    189                 void    GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
    190                 void    GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
    191                 void    GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    192                 void      GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
    193                 void      GetSolutionFromInputsEnthalpy(Vector<IssmDouble>* solutiong);
    194                 IssmDouble  GetStabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa);
     183                void             GetDofList(int** pdoflist,int approximation_enum,int setenum);
     184                void             GetVertexPidList(int* doflist);
     185                void           GetVertexSidList(int* sidlist);
     186                void           GetConnectivityList(int* connectivity);
     187                IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
     188                int            GetElementType(void);
     189                void           GetElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
     190                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
     191                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
     192                void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
     193                void             GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
     194                void             GetSolutionFromInputsEnthalpy(Vector<IssmDouble>* solutiong);
     195                IssmDouble     GetStabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa);
    195196                void    GetStrainRate3dPattyn(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
    196197                void    GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input);

  • issm/trunk-jpl/src/m/classes/groundingline.m

    r14361 r14605  
    4444                                        md = checkmessage(md,['bathymetry superior to bed on floating ice!']);
    4545                                end
    46                                 if strcmp(obj.migration,'SubelementMigration') & md.mesh.dimension==3,
    47                                         md = checkmessage(md,['SubelementMigration only implemented in 2d!']);
    48                                 end
    4946                        end
    5047

  • issm/trunk-jpl/src/m/classes/groundingline.py

    r14361 r14605  
    5959                        if any(md.geometry.bathymetry[pos]-md.geometry.bed[pos]>10**-9):
    6060                                md.checkmessage("bathymetry superior to bed on floating ice!")
    61                         if strcmp(self.migration,'SubelementMigration'):
    62                                 if md.mesh.dimension==3:
    63                                         md.checkmessage("SubelementMigration only implemented in 2d!")
    6461
    6562                return md
Note: See TracChangeset for help on using the changeset viewer.