Skip to content

Commit cfa6325

Browse files
committed
BUG:add 1e-10 to comparing thickness with crevasse to avoid machine epsilon
1 parent 9e86f27 commit cfa6325

File tree

2 files changed

+4
-4
lines changed

2 files changed

+4
-4
lines changed

src/c/analyses/LevelsetAnalysis.cpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -922,7 +922,7 @@ void LevelsetAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
922922
thickness_input->GetInputValue(&thickness,gauss);
923923
surface_input->GetInputValue(&surface,gauss);
924924

925-
if((surface_crevasse>surface || crevassedepth>crevasse_threshold*thickness) && bed<0.){
925+
if((surface_crevasse>surface || crevassedepth>(crevasse_threshold*thickness)-1e-10) && bed<0.){
926926
vec_constraint_nodes->SetValue(node->Pid(),1.0,INS_VAL);
927927
}
928928
}
@@ -972,7 +972,7 @@ void LevelsetAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
972972
surface_input->GetInputValue(&surface,gauss);
973973

974974
/*FIXME: not sure about levelset<0. && fabs(levelset)>-mig_max*dt! SHould maybe be distance<mig_max*dt*/
975-
if((surface_crevasse>surface || crevassedepth>crevasse_threshold*thickness) && bed<0. && levelset<0. && levelset>-mig_max*dt && constraint_nodes[node->Lid()]==0.){
975+
if((surface_crevasse>surface || crevassedepth>(crevasse_threshold*thickness-1e-10)) && bed<0. && levelset<0. && levelset>-mig_max*dt && constraint_nodes[node->Lid()]==0.){
976976
local_nflipped++;
977977
vec_constraint_nodes->SetValue(node->Pid(),1.0,INS_VAL);
978978
}

src/c/classes/Elements/Tria.cpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -562,13 +562,13 @@ void Tria::CalvingCrevasseDepth(){/*{{{*/
562562
strainrateeffective_input->GetInputValue(&straineffective,&gauss);
563563
n_input->GetInputValue(&n,&gauss);
564564
B_input->GetInputValue(&B,&gauss);
565-
if((straineffective <= 0.) || (thickness <= 0.) || (n<=0.)){
565+
if((straineffective <= 0.) || (thickness <= 0.) ){
566566
vH = 1e14;
567567
}
568568
else{
569569
vH = 0.5*B/thickness*pow(straineffective, (1./n)-1);
570570
}
571-
Kmax = 1 - 4.0*vH*(s1+s2+min(s1,s2))/(rho_ice*constant_g*(rho_seawater-rho_ice)/rho_seawater);
571+
Kmax = 1.0 - 4.0*vH*(s1+s2+min(s1,s2))/(rho_ice*constant_g*(rho_seawater-rho_ice)/rho_seawater);
572572
if(Kmax<0.) Kmax = 0.0;
573573
}
574574
else{

0 commit comments

Comments
 (0)