Changeset 26830
- Timestamp:
- 01/27/22 18:17:58 (3 years ago)
- 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 3416 3416 case StrainRateperpendicularEnum: this->StrainRateperpendicular(); break; 3417 3417 case SurfaceCrevasseEnum: this->CalvingCrevasseDepth(); break; 3418 case SigmaVMEnum: this->C alvingRateVonmises(); break;3418 case SigmaVMEnum: this->ComputeSigmaVM(); break; 3419 3419 case PartitioningEnum: this->inputs->SetInput(PartitioningEnum,this->lid,IssmComm::GetRank()); break; 3420 3420 } -
issm/trunk-jpl/src/c/classes/Elements/Element.h
r26810 r26830 238 238 virtual void ComputeDeviatoricStressTensor(void)=0; 239 239 virtual void ComputeSigmaNN(void)=0; 240 virtual void ComputeSigmaVM(void){_error_("not implemented yet");}; 240 241 virtual void ComputeStressTensor(void)=0; 241 242 virtual void ComputeEsaStrainAndVorticity(void)=0; -
issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
r26694 r26830 271 271 272 272 if(!this->IsOnBase()) return; 273 274 IssmDouble xyz_list[NUMVERTICES][3]; 275 IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ 273 this->ComputeSigmaVM(); 274 276 275 IssmDouble calvingratex[NUMVERTICES]; 277 276 IssmDouble calvingratey[NUMVERTICES]; 278 277 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; 286 281 287 282 /*Depth average B for stress calculation*/ 288 this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);289 283 this->InputDepthAverageAtBase(VxEnum,VxAverageEnum); 290 284 this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); … … 295 289 Input* gr_input = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input); 296 290 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);299 291 Input* smax_fl_input = this->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input); 300 292 Input* smax_gr_input = this->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input); 301 293 Input* sl_input = this->GetInput(SealevelEnum); _assert_(sl_input); 294 Input* sigma_vm_input = this->GetInput(SigmaVMEnum); _assert_(sigma_vm_input); 302 295 303 296 /* Start looping on the number of vertices: */ … … 307 300 308 301 /*Get velocity components and thickness*/ 309 B_input->GetInputValue(&B,&gauss);310 n_input->GetInputValue(&n,&gauss);311 302 vx_input->GetInputValue(&vx,&gauss); 312 303 vy_input->GetInputValue(&vy,&gauss); 304 sigma_vm_input->GetInputValue(&sigma_vm,&gauss); 313 305 gr_input->GetInputValue(&groundedice,&gauss); 314 306 bs_input->GetInputValue(&bed,&gauss); 315 307 smax_fl_input->GetInputValue(&sigma_max_floating,&gauss); 316 308 smax_gr_input->GetInputValue(&sigma_max_grounded,&gauss); 317 vel=sqrt(vx*vx+vy*vy)+1.e-14;318 309 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));335 310 336 311 /*Tensile stress threshold*/ … … 346 321 } 347 322 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; 350 325 } 351 326 calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]); … … 356 331 this->AddBasalInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum); 357 332 this->AddBasalInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); 358 this->AddBasalInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum);359 333 360 334 this->InputExtrude(CalvingratexEnum,-1); 361 335 this->InputExtrude(CalvingrateyEnum,-1); 362 336 this->InputExtrude(CalvingCalvingrateEnum,-1); 363 this->InputExtrude(SigmaVMEnum,-1);364 337 } 365 338 /*}}}*/ … … 847 820 this->AddInput(StressTensoryzEnum,&sigma_yz[0],P1DGEnum); 848 821 this->AddInput(StressTensorzzEnum,&sigma_zz[0],P1DGEnum); 822 } 823 /*}}}*/ 824 void 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); 849 880 } 850 881 /*}}}*/ -
issm/trunk-jpl/src/c/classes/Elements/Penta.h
r26800 r26830 49 49 void AddControlInput(int input_enum,Inputs* inputs,IoModel* iomodel,IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max, int interpolation_enum,int id); 50 50 void ControlInputExtrude(int enum_type,int start); 51 void ComputeSigmaVM(void); 51 52 void DatasetInputExtrude(int enum_type,int start); 52 53 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 313 313 void Tria::CalvingRateVonmises(){/*{{{*/ 314 314 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*/ 317 319 IssmDouble calvingratex[NUMVERTICES]; 318 320 IssmDouble calvingratey[NUMVERTICES]; 319 321 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; 327 325 328 326 /*Retrieve all inputs and parameters we will need*/ 329 327 Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input); 330 328 Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input); 331 Input* B_input = this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input);332 329 Input* gr_input = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input); 333 330 Input* bs_input = this->GetInput(BaseEnum); _assert_(bs_input); 334 331 Input* smax_fl_input = this->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input); 335 332 Input* smax_gr_input = this->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input); 336 Input* n_input = this->GetInput(MaterialsRheologyNEnum); _assert_(n_input);337 333 Input* sl_input = this->GetInput(SealevelEnum); _assert_(sl_input); 334 Input* sigma_vm_input = this->GetInput(SigmaVMEnum); _assert_(sigma_vm_input); 338 335 339 336 /* Start looping on the number of vertices: */ … … 343 340 344 341 /*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); 347 343 vx_input->GetInputValue(&vx,&gauss); 348 344 vy_input->GetInputValue(&vy,&gauss); … … 351 347 smax_fl_input->GetInputValue(&sigma_max_floating,&gauss); 352 348 smax_gr_input->GetInputValue(&sigma_max_grounded,&gauss); 353 vel=sqrt(vx*vx+vy*vy)+1.e-14;354 349 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 test374 //sigma_max = 1000.e+3; //GRL375 //if(groundedice<0) sigma_max=150.e+3;376 350 377 351 /*Tensile stress threshold*/ … … 381 355 sigma_max = sigma_max_grounded; 382 356 /*Assign values*/ 383 if(bed>sealevel){ // Changed 0. to sealevel357 if(bed>sealevel){ 384 358 calvingratex[iv]=0.; 385 359 calvingratey[iv]=0.; 386 360 } 387 361 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; 390 364 } 391 365 calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]); 392 393 366 } 394 367 … … 397 370 this->AddInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum); 398 371 this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); 399 this->AddInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum);400 372 } 401 373 /*}}}*/ … … 947 919 } 948 920 /*}}}*/ 949 void Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/921 void Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/ 950 922 951 923 IssmDouble xyz_list[NUMVERTICES][3]; … … 1043 1015 delete gauss; 1044 1016 } 1017 } 1018 /*}}}*/ 1019 void 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); 1045 1067 } 1046 1068 /*}}}*/ -
issm/trunk-jpl/src/c/classes/Elements/Tria.h
r26800 r26830 65 65 void ComputeEsaStrainAndVorticity(); 66 66 void ComputeSigmaNN(); 67 void ComputeSigmaVM(); 67 68 void ComputeStressTensor(); 68 69 void ComputeSurfaceNormalVelocity();
Note:
See TracChangeset
for help on using the changeset viewer.