Changeset 26830


Ignore:
Timestamp:
01/27/22 18:17:58 (3 years ago)
Author:
Mathieu Morlighem
Message:

CHG: added sigmavm for Julia

Location:
issm/trunk-jpl/src/c/classes/Elements
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r26814 r26830  
    34163416                case StrainRateperpendicularEnum: this->StrainRateperpendicular(); break;
    34173417                case SurfaceCrevasseEnum: this->CalvingCrevasseDepth(); break;
    3418                 case SigmaVMEnum: this->CalvingRateVonmises(); break;
     3418                case SigmaVMEnum: this->ComputeSigmaVM(); break;
    34193419                case PartitioningEnum: this->inputs->SetInput(PartitioningEnum,this->lid,IssmComm::GetRank()); break;
    34203420        }

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

    r26810 r26830  
    238238                virtual void       ComputeDeviatoricStressTensor(void)=0;
    239239                virtual void       ComputeSigmaNN(void)=0;
     240                virtual void       ComputeSigmaVM(void){_error_("not implemented yet");};
    240241                virtual void       ComputeStressTensor(void)=0;
    241242                virtual void       ComputeEsaStrainAndVorticity(void)=0;

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

    r26694 r26830  
    271271
    272272        if(!this->IsOnBase()) return;
    273 
    274         IssmDouble  xyz_list[NUMVERTICES][3];
    275         IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     273        this->ComputeSigmaVM();
     274
    276275        IssmDouble  calvingratex[NUMVERTICES];
    277276        IssmDouble  calvingratey[NUMVERTICES];
    278277        IssmDouble  calvingrate[NUMVERTICES];
    279         IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
    280         IssmDouble  B,sigma_max,sigma_max_floating,sigma_max_grounded,n;
    281         IssmDouble  epse_2,groundedice,bed,sealevel;
    282         IssmDouble  sigma_vm[NUMVERTICES];
    283 
    284         /* Get node coordinates and dof list: */
    285         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     278        IssmDouble  vx,vy;
     279        IssmDouble  sigma_vm,sigma_max,sigma_max_floating,sigma_max_grounded;
     280        IssmDouble  groundedice,bed,sealevel;
    286281
    287282        /*Depth average B for stress calculation*/
    288         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
    289283        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
    290284        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
    … …  
    295289        Input* gr_input = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
    296290        Input* bs_input = this->GetInput(BaseEnum);                    _assert_(bs_input);
    297         Input* B_input  = this->GetInput(MaterialsRheologyBbarEnum);   _assert_(B_input);
    298         Input* n_input  = this->GetInput(MaterialsRheologyNEnum);   _assert_(n_input);
    299291        Input* smax_fl_input = this->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input);
    300292        Input* smax_gr_input = this->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input);
    301293        Input* sl_input  = this->GetInput(SealevelEnum); _assert_(sl_input);
     294        Input* sigma_vm_input = this->GetInput(SigmaVMEnum); _assert_(sigma_vm_input);
    302295
    303296        /* Start looping on the number of vertices: */
    … …  
    307300
    308301                /*Get velocity components and thickness*/
    309                 B_input->GetInputValue(&B,&gauss);
    310                 n_input->GetInputValue(&n,&gauss);
    311302                vx_input->GetInputValue(&vx,&gauss);
    312303                vy_input->GetInputValue(&vy,&gauss);
     304                sigma_vm_input->GetInputValue(&sigma_vm,&gauss);
    313305                gr_input->GetInputValue(&groundedice,&gauss);
    314306                bs_input->GetInputValue(&bed,&gauss);
    315307                smax_fl_input->GetInputValue(&sigma_max_floating,&gauss);
    316308                smax_gr_input->GetInputValue(&sigma_max_grounded,&gauss);
    317                 vel=sqrt(vx*vx+vy*vy)+1.e-14;
    318309                sl_input->GetInputValue(&sealevel,&gauss);
    319 
    320                 /*Compute strain rate and viscosity: */
    321                 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input);
    322 
    323                 /*Get Eigen values*/
    324                 Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]);
    325                 _assert_(!xIsNan<IssmDouble>(lambda1));
    326                 _assert_(!xIsNan<IssmDouble>(lambda2));
    327 
    328                 /*Process Eigen values (only account for extension)*/
    329                 lambda1 = max(lambda1,0.);
    330                 lambda2 = max(lambda2,0.);
    331 
    332                 /*Calculate sigma_vm*/
    333                 epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
    334                 sigma_vm[iv] = sqrt(3.) * B * pow(epse_2,1./(2.*n));
    335310
    336311                /*Tensile stress threshold*/
    … …  
    346321                }
    347322                else{
    348                         calvingratex[iv]=vx*sigma_vm[iv]/sigma_max;
    349                         calvingratey[iv]=vy*sigma_vm[iv]/sigma_max;
     323                        calvingratex[iv]=vx*sigma_vm/sigma_max;
     324                        calvingratey[iv]=vy*sigma_vm/sigma_max;
    350325                }
    351326                calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
    … …  
    356331        this->AddBasalInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
    357332        this->AddBasalInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
    358         this->AddBasalInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum);
    359333
    360334        this->InputExtrude(CalvingratexEnum,-1);
    361335        this->InputExtrude(CalvingrateyEnum,-1);
    362336        this->InputExtrude(CalvingCalvingrateEnum,-1);
    363         this->InputExtrude(SigmaVMEnum,-1);
    364337}
    365338/*}}}*/
    … …  
    847820        this->AddInput(StressTensoryzEnum,&sigma_yz[0],P1DGEnum);
    848821        this->AddInput(StressTensorzzEnum,&sigma_zz[0],P1DGEnum);
     822}
     823/*}}}*/
     824void       Penta::ComputeSigmaVM(){/*{{{*/
     825
     826        if(!this->IsOnBase()) return;
     827
     828        IssmDouble  xyz_list[NUMVERTICES][3];
     829        IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     830        IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
     831        IssmDouble  B,n;
     832        IssmDouble  sigma_vm[NUMVERTICES];
     833
     834        /* Get node coordinates and dof list: */
     835        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     836
     837        /*Depth average B for stress calculation*/
     838        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
     839        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
     840        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
     841
     842        /*Retrieve all inputs and parameters we will need*/
     843        Input* vx_input = this->GetInput(VxAverageEnum); _assert_(vx_input);
     844        Input* vy_input = this->GetInput(VyAverageEnum); _assert_(vy_input);
     845        Input* B_input  = this->GetInput(MaterialsRheologyBbarEnum);   _assert_(B_input);
     846        Input* n_input  = this->GetInput(MaterialsRheologyNEnum);   _assert_(n_input);
     847
     848        /* Start looping on the number of vertices: */
     849        GaussPenta gauss;
     850        for(int iv=0;iv<3;iv++){
     851                gauss.GaussVertex(iv);
     852
     853                /*Get velocity components and thickness*/
     854                B_input->GetInputValue(&B,&gauss);
     855                n_input->GetInputValue(&n,&gauss);
     856                vx_input->GetInputValue(&vx,&gauss);
     857                vy_input->GetInputValue(&vy,&gauss);
     858                vel=sqrt(vx*vx+vy*vy)+1.e-14;
     859
     860                /*Compute strain rate and viscosity: */
     861                this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input);
     862
     863                /*Get Eigen values*/
     864                Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]);
     865                _assert_(!xIsNan<IssmDouble>(lambda1));
     866                _assert_(!xIsNan<IssmDouble>(lambda2));
     867
     868                /*Process Eigen values (only account for extension)*/
     869                lambda1 = max(lambda1,0.);
     870                lambda2 = max(lambda2,0.);
     871
     872                /*Calculate sigma_vm*/
     873                IssmDouble epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
     874                sigma_vm[iv] = sqrt(3.) * B * pow(epse_2,1./(2.*n));
     875        }
     876
     877        /*Add input*/
     878        this->AddBasalInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum);
     879        this->InputExtrude(SigmaVMEnum,-1);
    849880}
    850881/*}}}*/

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

    r26800 r26830  
    4949                void           AddControlInput(int input_enum,Inputs* inputs,IoModel* iomodel,IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max, int interpolation_enum,int id);
    5050                void           ControlInputExtrude(int enum_type,int start);
     51                void           ComputeSigmaVM(void);
    5152                void           DatasetInputExtrude(int enum_type,int start);
    5253                void           DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,Inputs* inputs,IoModel* iomodel,int input_enum);

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

    r26800 r26830  
    313313void       Tria::CalvingRateVonmises(){/*{{{*/
    314314
    315         IssmDouble  xyz_list[NUMVERTICES][3];
    316         IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     315        /*First, compute Von Mises Stress*/
     316        this->ComputeSigmaVM();
     317
     318        /*Now compute calving rate*/
    317319        IssmDouble  calvingratex[NUMVERTICES];
    318320        IssmDouble  calvingratey[NUMVERTICES];
    319321        IssmDouble  calvingrate[NUMVERTICES];
    320         IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
    321         IssmDouble  sigma_vm[NUMVERTICES];
    322         IssmDouble  B,sigma_max,sigma_max_floating,sigma_max_grounded,n;
    323         IssmDouble  epse_2,groundedice,bed,sealevel;            // added sealevel
    324 
    325         /* Get node coordinates and dof list: */
    326         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     322        IssmDouble  sigma_vm,vx,vy;
     323        IssmDouble  sigma_max,sigma_max_floating,sigma_max_grounded,n;
     324        IssmDouble  groundedice,bed,sealevel;
    327325
    328326        /*Retrieve all inputs and parameters we will need*/
    329327        Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input);
    330328        Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input);
    331         Input* B_input  = this->GetInput(MaterialsRheologyBbarEnum);   _assert_(B_input);
    332329        Input* gr_input = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
    333330        Input* bs_input = this->GetInput(BaseEnum);                    _assert_(bs_input);
    334331        Input* smax_fl_input = this->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input);
    335332        Input* smax_gr_input = this->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input);
    336         Input* n_input  = this->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
    337333        Input* sl_input  = this->GetInput(SealevelEnum); _assert_(sl_input);
     334        Input* sigma_vm_input  = this->GetInput(SigmaVMEnum); _assert_(sigma_vm_input);
    338335
    339336        /* Start looping on the number of vertices: */
    … …  
    343340
    344341                /*Get velocity components and thickness*/
    345                 B_input->GetInputValue(&B,&gauss);
    346                 n_input->GetInputValue(&n,&gauss);
     342                sigma_vm_input->GetInputValue(&sigma_vm,&gauss);
    347343                vx_input->GetInputValue(&vx,&gauss);
    348344                vy_input->GetInputValue(&vy,&gauss);
    … …  
    351347                smax_fl_input->GetInputValue(&sigma_max_floating,&gauss);
    352348                smax_gr_input->GetInputValue(&sigma_max_grounded,&gauss);
    353                 vel=sqrt(vx*vx+vy*vy)+1.e-14;
    354349                sl_input->GetInputValue(&sealevel,&gauss);
    355 
    356                 /*Compute strain rate and viscosity: */
    357                 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input);
    358 
    359                 /*Get Eigen values*/
    360                 Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]);
    361                 _assert_(!xIsNan<IssmDouble>(lambda1));
    362                 _assert_(!xIsNan<IssmDouble>(lambda2));
    363 
    364                 /*Process Eigen values (only account for extension)*/
    365                 lambda1 = max(lambda1,0.);
    366                 lambda2 = max(lambda2,0.);
    367 
    368                 /*Calculate sigma_vm*/
    369                 epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
    370                 sigma_vm[iv]  = sqrt(3.) * B * pow(epse_2,1./(2.*n));
    371 
    372                 /*OLD (keep for a little bit)*/
    373                 //sigma_max = 800.e+3; //IUGG previous test
    374                 //sigma_max = 1000.e+3; //GRL
    375                 //if(groundedice<0) sigma_max=150.e+3;
    376350
    377351                /*Tensile stress threshold*/
    … …  
    381355                 sigma_max = sigma_max_grounded;
    382356                /*Assign values*/
    383                 if(bed>sealevel){               // Changed 0. to sealevel
     357                if(bed>sealevel){
    384358                        calvingratex[iv]=0.;
    385359                        calvingratey[iv]=0.;
    386360                }
    387361                else{
    388                         calvingratex[iv]=vx*sigma_vm[iv]/sigma_max;
    389                         calvingratey[iv]=vy*sigma_vm[iv]/sigma_max;
     362                        calvingratex[iv]=vx*sigma_vm/sigma_max;
     363                        calvingratey[iv]=vy*sigma_vm/sigma_max;
    390364                }
    391365                calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
    392 
    393366        }
    394367
    … …  
    397370        this->AddInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
    398371        this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
    399         this->AddInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum);
    400372}
    401373/*}}}*/
    … …  
    947919}
    948920/*}}}*/
    949 void         Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/
     921void          Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/
    950922
    951923        IssmDouble  xyz_list[NUMVERTICES][3];
    … …  
    10431015                delete gauss;
    10441016        }
     1017}
     1018/*}}}*/
     1019void       Tria::ComputeSigmaVM(){/*{{{*/
     1020
     1021        IssmDouble  xyz_list[NUMVERTICES][3];
     1022        IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     1023        IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
     1024        IssmDouble  sigma_vm[NUMVERTICES];
     1025        IssmDouble  B,n;
     1026
     1027        /* Get node coordinates and dof list: */
     1028        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     1029
     1030        /*Retrieve all inputs and parameters we will need*/
     1031        Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input);
     1032        Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input);
     1033        Input* B_input  = this->GetInput(MaterialsRheologyBbarEnum);   _assert_(B_input);
     1034        Input* n_input  = this->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
     1035
     1036        /* Start looping on the number of vertices: */
     1037        GaussTria gauss;
     1038        for(int iv=0;iv<NUMVERTICES;iv++){
     1039                gauss.GaussVertex(iv);
     1040
     1041                /*Get velocity components and thickness*/
     1042                B_input->GetInputValue(&B,&gauss);
     1043                n_input->GetInputValue(&n,&gauss);
     1044                vx_input->GetInputValue(&vx,&gauss);
     1045                vy_input->GetInputValue(&vy,&gauss);
     1046                vel=sqrt(vx*vx+vy*vy)+1.e-14;
     1047
     1048                /*Compute strain rate and viscosity: */
     1049                this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input);
     1050
     1051                /*Get Eigen values*/
     1052                Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]);
     1053                _assert_(!xIsNan<IssmDouble>(lambda1));
     1054                _assert_(!xIsNan<IssmDouble>(lambda2));
     1055
     1056                /*Process Eigen values (only account for extension)*/
     1057                lambda1 = max(lambda1,0.);
     1058                lambda2 = max(lambda2,0.);
     1059
     1060                /*Calculate sigma_vm*/
     1061                IssmDouble epse_2 = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
     1062                sigma_vm[iv]  = sqrt(3.) * B * pow(epse_2,1./(2.*n));
     1063        }
     1064
     1065        /*Add input*/
     1066        this->AddInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum);
    10451067}
    10461068/*}}}*/

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

    r26800 r26830  
    6565                void        ComputeEsaStrainAndVorticity();
    6666                void        ComputeSigmaNN();
     67                void        ComputeSigmaVM();
    6768                void        ComputeStressTensor();
    6869                void        ComputeSurfaceNormalVelocity();
Note: See TracChangeset for help on using the changeset viewer.