Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26829) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26830) @@ -3416,5 +3416,5 @@ case StrainRateperpendicularEnum: this->StrainRateperpendicular(); break; case SurfaceCrevasseEnum: this->CalvingCrevasseDepth(); break; - case SigmaVMEnum: this->CalvingRateVonmises(); break; + case SigmaVMEnum: this->ComputeSigmaVM(); break; case PartitioningEnum: this->inputs->SetInput(PartitioningEnum,this->lid,IssmComm::GetRank()); break; } Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26829) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26830) @@ -238,4 +238,5 @@ virtual void ComputeDeviatoricStressTensor(void)=0; virtual void ComputeSigmaNN(void)=0; + virtual void ComputeSigmaVM(void){_error_("not implemented yet");}; virtual void ComputeStressTensor(void)=0; virtual void ComputeEsaStrainAndVorticity(void)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 26829) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 26830) @@ -271,20 +271,14 @@ if(!this->IsOnBase()) return; - - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + this->ComputeSigmaVM(); + IssmDouble calvingratex[NUMVERTICES]; IssmDouble calvingratey[NUMVERTICES]; IssmDouble calvingrate[NUMVERTICES]; - IssmDouble lambda1,lambda2,ex,ey,vx,vy,vel; - IssmDouble B,sigma_max,sigma_max_floating,sigma_max_grounded,n; - IssmDouble epse_2,groundedice,bed,sealevel; - IssmDouble sigma_vm[NUMVERTICES]; - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + IssmDouble vx,vy; + IssmDouble sigma_vm,sigma_max,sigma_max_floating,sigma_max_grounded; + IssmDouble groundedice,bed,sealevel; /*Depth average B for stress calculation*/ - this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); this->InputDepthAverageAtBase(VxEnum,VxAverageEnum); this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); @@ -295,9 +289,8 @@ Input* gr_input = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input); Input* bs_input = this->GetInput(BaseEnum); _assert_(bs_input); - Input* B_input = this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input); - Input* n_input = this->GetInput(MaterialsRheologyNEnum); _assert_(n_input); Input* smax_fl_input = this->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input); Input* smax_gr_input = this->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input); Input* sl_input = this->GetInput(SealevelEnum); _assert_(sl_input); + Input* sigma_vm_input = this->GetInput(SigmaVMEnum); _assert_(sigma_vm_input); /* Start looping on the number of vertices: */ @@ -307,30 +300,12 @@ /*Get velocity components and thickness*/ - B_input->GetInputValue(&B,&gauss); - n_input->GetInputValue(&n,&gauss); vx_input->GetInputValue(&vx,&gauss); vy_input->GetInputValue(&vy,&gauss); + sigma_vm_input->GetInputValue(&sigma_vm,&gauss); gr_input->GetInputValue(&groundedice,&gauss); bs_input->GetInputValue(&bed,&gauss); smax_fl_input->GetInputValue(&sigma_max_floating,&gauss); smax_gr_input->GetInputValue(&sigma_max_grounded,&gauss); - vel=sqrt(vx*vx+vy*vy)+1.e-14; sl_input->GetInputValue(&sealevel,&gauss); - - /*Compute strain rate and viscosity: */ - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); - - /*Get Eigen values*/ - Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]); - _assert_(!xIsNan(lambda1)); - _assert_(!xIsNan(lambda2)); - - /*Process Eigen values (only account for extension)*/ - lambda1 = max(lambda1,0.); - lambda2 = max(lambda2,0.); - - /*Calculate sigma_vm*/ - epse_2 = 1./2. *(lambda1*lambda1 + lambda2*lambda2); - sigma_vm[iv] = sqrt(3.) * B * pow(epse_2,1./(2.*n)); /*Tensile stress threshold*/ @@ -346,6 +321,6 @@ } else{ - calvingratex[iv]=vx*sigma_vm[iv]/sigma_max; - calvingratey[iv]=vy*sigma_vm[iv]/sigma_max; + calvingratex[iv]=vx*sigma_vm/sigma_max; + calvingratey[iv]=vy*sigma_vm/sigma_max; } calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]); @@ -356,10 +331,8 @@ this->AddBasalInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum); this->AddBasalInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); - this->AddBasalInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum); this->InputExtrude(CalvingratexEnum,-1); this->InputExtrude(CalvingrateyEnum,-1); this->InputExtrude(CalvingCalvingrateEnum,-1); - this->InputExtrude(SigmaVMEnum,-1); } /*}}}*/ @@ -847,4 +820,62 @@ this->AddInput(StressTensoryzEnum,&sigma_yz[0],P1DGEnum); this->AddInput(StressTensorzzEnum,&sigma_zz[0],P1DGEnum); +} +/*}}}*/ +void Penta::ComputeSigmaVM(){/*{{{*/ + + if(!this->IsOnBase()) return; + + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble lambda1,lambda2,ex,ey,vx,vy,vel; + IssmDouble B,n; + IssmDouble sigma_vm[NUMVERTICES]; + + /* Get node coordinates and dof list: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*Depth average B for stress calculation*/ + this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); + this->InputDepthAverageAtBase(VxEnum,VxAverageEnum); + this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); + + /*Retrieve all inputs and parameters we will need*/ + Input* vx_input = this->GetInput(VxAverageEnum); _assert_(vx_input); + Input* vy_input = this->GetInput(VyAverageEnum); _assert_(vy_input); + Input* B_input = this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input); + Input* n_input = this->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + + /* Start looping on the number of vertices: */ + GaussPenta gauss; + for(int iv=0;iv<3;iv++){ + gauss.GaussVertex(iv); + + /*Get velocity components and thickness*/ + B_input->GetInputValue(&B,&gauss); + n_input->GetInputValue(&n,&gauss); + vx_input->GetInputValue(&vx,&gauss); + vy_input->GetInputValue(&vy,&gauss); + vel=sqrt(vx*vx+vy*vy)+1.e-14; + + /*Compute strain rate and viscosity: */ + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); + + /*Get Eigen values*/ + Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]); + _assert_(!xIsNan(lambda1)); + _assert_(!xIsNan(lambda2)); + + /*Process Eigen values (only account for extension)*/ + lambda1 = max(lambda1,0.); + lambda2 = max(lambda2,0.); + + /*Calculate sigma_vm*/ + IssmDouble epse_2 = 1./2. *(lambda1*lambda1 + lambda2*lambda2); + sigma_vm[iv] = sqrt(3.) * B * pow(epse_2,1./(2.*n)); + } + + /*Add input*/ + this->AddBasalInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum); + this->InputExtrude(SigmaVMEnum,-1); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26829) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26830) @@ -49,4 +49,5 @@ void AddControlInput(int input_enum,Inputs* inputs,IoModel* iomodel,IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max, int interpolation_enum,int id); void ControlInputExtrude(int enum_type,int start); + void ComputeSigmaVM(void); void DatasetInputExtrude(int enum_type,int start); void DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,Inputs* inputs,IoModel* iomodel,int input_enum); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26829) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26830) @@ -313,27 +313,24 @@ void Tria::CalvingRateVonmises(){/*{{{*/ - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + /*First, compute Von Mises Stress*/ + this->ComputeSigmaVM(); + + /*Now compute calving rate*/ IssmDouble calvingratex[NUMVERTICES]; IssmDouble calvingratey[NUMVERTICES]; IssmDouble calvingrate[NUMVERTICES]; - IssmDouble lambda1,lambda2,ex,ey,vx,vy,vel; - IssmDouble sigma_vm[NUMVERTICES]; - IssmDouble B,sigma_max,sigma_max_floating,sigma_max_grounded,n; - IssmDouble epse_2,groundedice,bed,sealevel; // added sealevel - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + IssmDouble sigma_vm,vx,vy; + IssmDouble sigma_max,sigma_max_floating,sigma_max_grounded,n; + IssmDouble groundedice,bed,sealevel; /*Retrieve all inputs and parameters we will need*/ Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input); Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input); - Input* B_input = this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input); Input* gr_input = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input); Input* bs_input = this->GetInput(BaseEnum); _assert_(bs_input); Input* smax_fl_input = this->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input); Input* smax_gr_input = this->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input); - Input* n_input = this->GetInput(MaterialsRheologyNEnum); _assert_(n_input); Input* sl_input = this->GetInput(SealevelEnum); _assert_(sl_input); + Input* sigma_vm_input = this->GetInput(SigmaVMEnum); _assert_(sigma_vm_input); /* Start looping on the number of vertices: */ @@ -343,6 +340,5 @@ /*Get velocity components and thickness*/ - B_input->GetInputValue(&B,&gauss); - n_input->GetInputValue(&n,&gauss); + sigma_vm_input->GetInputValue(&sigma_vm,&gauss); vx_input->GetInputValue(&vx,&gauss); vy_input->GetInputValue(&vy,&gauss); @@ -351,27 +347,5 @@ smax_fl_input->GetInputValue(&sigma_max_floating,&gauss); smax_gr_input->GetInputValue(&sigma_max_grounded,&gauss); - vel=sqrt(vx*vx+vy*vy)+1.e-14; sl_input->GetInputValue(&sealevel,&gauss); - - /*Compute strain rate and viscosity: */ - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); - - /*Get Eigen values*/ - Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]); - _assert_(!xIsNan(lambda1)); - _assert_(!xIsNan(lambda2)); - - /*Process Eigen values (only account for extension)*/ - lambda1 = max(lambda1,0.); - lambda2 = max(lambda2,0.); - - /*Calculate sigma_vm*/ - epse_2 = 1./2. *(lambda1*lambda1 + lambda2*lambda2); - sigma_vm[iv] = sqrt(3.) * B * pow(epse_2,1./(2.*n)); - - /*OLD (keep for a little bit)*/ - //sigma_max = 800.e+3; //IUGG previous test - //sigma_max = 1000.e+3; //GRL - //if(groundedice<0) sigma_max=150.e+3; /*Tensile stress threshold*/ @@ -381,14 +355,13 @@ sigma_max = sigma_max_grounded; /*Assign values*/ - if(bed>sealevel){ // Changed 0. to sealevel + if(bed>sealevel){ calvingratex[iv]=0.; calvingratey[iv]=0.; } else{ - calvingratex[iv]=vx*sigma_vm[iv]/sigma_max; - calvingratey[iv]=vy*sigma_vm[iv]/sigma_max; + calvingratex[iv]=vx*sigma_vm/sigma_max; + calvingratey[iv]=vy*sigma_vm/sigma_max; } calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]); - } @@ -397,5 +370,4 @@ this->AddInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum); this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); - this->AddInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum); } /*}}}*/ @@ -947,5 +919,5 @@ } /*}}}*/ -void Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/ +void Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/ IssmDouble xyz_list[NUMVERTICES][3]; @@ -1043,4 +1015,54 @@ delete gauss; } +} +/*}}}*/ +void Tria::ComputeSigmaVM(){/*{{{*/ + + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble lambda1,lambda2,ex,ey,vx,vy,vel; + IssmDouble sigma_vm[NUMVERTICES]; + IssmDouble B,n; + + /* Get node coordinates and dof list: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*Retrieve all inputs and parameters we will need*/ + Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input); + Input* B_input = this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input); + Input* n_input = this->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + + /* Start looping on the number of vertices: */ + GaussTria gauss; + for(int iv=0;ivGetInputValue(&B,&gauss); + n_input->GetInputValue(&n,&gauss); + vx_input->GetInputValue(&vx,&gauss); + vy_input->GetInputValue(&vy,&gauss); + vel=sqrt(vx*vx+vy*vy)+1.e-14; + + /*Compute strain rate and viscosity: */ + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); + + /*Get Eigen values*/ + Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]); + _assert_(!xIsNan(lambda1)); + _assert_(!xIsNan(lambda2)); + + /*Process Eigen values (only account for extension)*/ + lambda1 = max(lambda1,0.); + lambda2 = max(lambda2,0.); + + /*Calculate sigma_vm*/ + IssmDouble epse_2 = 1./2. *(lambda1*lambda1 + lambda2*lambda2); + sigma_vm[iv] = sqrt(3.) * B * pow(epse_2,1./(2.*n)); + } + + /*Add input*/ + this->AddInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26829) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26830) @@ -65,4 +65,5 @@ void ComputeEsaStrainAndVorticity(); void ComputeSigmaNN(); + void ComputeSigmaVM(); void ComputeStressTensor(); void ComputeSurfaceNormalVelocity();